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o  describe  the  dynamic  evolution  of  frazil  ice  in  turbulent  natural  water  bodies,  the  basic  equation  for  dynamic  frazil 
crystal  number  continuity  and  the  basic  equation  of  heat  balance  for  a  differential  volume  are  developed.  Crystal  growth 
and  nucleation  of  new  crystals  are  the  major  parameters  in  these  equations.  Expressions  tor  the  growth  rate  ol  frazil  ice 
crystals  are  described.  The  growth  rate  along  the  major  axis  is  controlled  by  heat  transfer.  The  heat  transfer  coefficient 
is  a  function  of  crystal  size,  the  fluid  turbulence,  and  the  fluid  properties.  The  magnitude  of  inertial  and  buoyancy  forces 
on  the  ice  crystals  are  determined  as  is  their  influence  on  the  heat  transfer.  Spontaneous  nucleation  of  ice  can  be  discounter 
secondary  nucleation  is  responsible  for  the  vast  majority  of  frazil  ice  crystals.  The  theoretical  rate  of  secondary  nuclea¬ 
tion  is  partially  modeled  as  a  function  of  the  supercooling,  fluid  turbulence  and  crystal  size  distribution.  A  simple  an- 
alytical  solution  of  the  basic  equations  is  developed  for  the  growth  of  frazil  ice  in  a  well-mixed,  steady-state  crystallizer. 
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This  monograph  provides  a  technical  analysis  of  the  state  of  thv  rt  in  some  facets  of  the  phys¬ 
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Civil  Engineering,  R.M.  Parsons  Laboratory  for  Water  Resources  and  Hydrodynamics.  Massachu¬ 
setts  Institute  of  Technology .  It  was  prepared  by  the  author  at  the  Massachusetts  Institute  of  Tech 
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fillment  of  the  requirements  for  the  Degree  of  Master  of  Science  in  Civil  Engineering.  Funding  was 
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commercial  products. 
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NOMENCLATURE 


a  Constant 

a1.a3.a3  Constants 

a',  a"  Linearized  intrinsic  kinetic  growth  rate  constants 
Acceleration  of  fluid 

<  Relative  acceleration  of  crystal  and  fluid 

Surface  area  of  disk 

r  Total  surface  area  of  crystals  in  R 

r)  Area  of  crystal  of  size  r 

Constant 
Birth  function 
Heat  capacity  of  fluid 
,  Drag  coefficient 

r  Fluid  impurity  concentration 

Heat  capacity  of  ice 
Death  function 

b  Collision  energy  created  by  collisions  of  crystals  and  boundaries 

Pure  straining  motion  of  fluid 
£w  Fluid  strain  rate  along  crystal  rotation  axis  of  symmetry 

F(r  1 .  r2 )  Collision  energy  created  by  collisions  of  crystals  of  size  rt  and  r2 

Et  Rate  of  energy  transfer  by  collision 

F,  Number  of  particles  generated  per  unit  of  collision  energy 

Fj  Fraction  of  particles  surviving  to  become  crystals 

F0  Drag  force 

g  Gravity 

g  Reduced  gravity 

G  Crystal  growth  rate 

Gk  Kinetic  growth  rate 

h  Heat  transfer  coefficient 

he  Heat  transfer  coefficient  from  disk  edge 

h{  Heat  transfer  coefficient  from  disk  face 

k  Thermal  conductivity  of  water 

A.'a  Geometric  shape  factor 

A'v  Volumetric  shape  factor 

/.  Mean  latent  heat  of  fusion  of  ice 

l.te  Latent  heat  of  fusion  of  ice  as  a  function  of  the  equilibrium  temperature 

m*  Nondimensional  crystal  size 

m(r)  Mass  of  crystal  of  size  r 

Mj  Total  mass  of  crystals  in  R 

n  Size  distribution  function 
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FRAZIL  ICE  DYNAMICS 


Steven  F.  Daly 


INTRODUCTION 

Frazil  ice,  fine  spicule,  plate  or  discoid  crystals  in  supercooled  turbulent  water,  is  commonplace 
on  many  northern  rivers  and  lakes  during  the  winter,  but  the  processes  that  create  it  are  not  well 
understood.  Frazil  is  formed  when  heat  is  withdrawn  from  a  body  of  turbulent  water  that  is  at  the 
freezing  point.  The  temperature  of  the  water  follows  a  qualitatively  well-known  sequence:  it  falls 
below  the  freezing  temperature  to  a  minimum  and  then  returns  to  the  freezing  point.  This  temper¬ 
ature  sequence  represents  a  dynamic  balance  between  the  latent  heat  released  by  the  growing  frazil 
ice  crystals  and  the  heat  lost  to  the  environment  from  the  water. 

Frazil  ice  can  block  water  supply  intakes,  hydroelectric  plant  intakes,  irrigation  and  water  supply 
canals,  and  it  can  form  ice  jams  that  can  hlock  an  entire  river  cross  section  and  cause  extensive 
Hooding.  Although  the  adverse  effects  of  frazil  cost  the  world  uncounted  millions  of  dollars  each 
year,  currently  no  quantitative  estimates  about  any  aspect  of  this  sequence  of  frazil  formation  can 
be  made.  The  deficiencies  in  the  current  state  of  knowledge  are  “severely  hindering  development 
of  rational  design  methods  for  avoidance  or  alleviation  of  frazil  ice  problems"  (ASCF  Task  Com¬ 
mittee  1974). 

There  exists  then  an  obvious  need  for  a  comprehensive,  physically  based  quantitative  model  of 
the  basic  processes  that  control  and  direct  the  growth  of  frazil  in  natural  water  bodies.  The  slow 
development  of  such  a  model  has  been  frustrating  to  the  engineers  who  must  cope  with  the  unique 
problems  of  frazil.  Despite  such  provocative  statements  as  “the  phenomenon  of  frazil...  is  similar 
to  the  more  general  one  of  crystallization  in  a  supersaturated  medium"  (Michel  1963),  very  seldom 
has  an  attempt  been  made  to  study  frazil  in  the  wider  context  of  crystal  growth  from  a  solution. 
The  problems  of  large-scale  industrial  crystallization,  however,  can  provide  insight  into  the  frazil 
problem. 

The  goal  of  this  report  is  to  take  the  first  steps  toward  development  of  a  comprehensive,  quanti¬ 
tative  model  of  the  process  of  frazil  formation  by  relating  knowledge  of  industrial  crystallization  to 
the  conditions  present  in  natural  water  bodies.  Toward  this  end.  this  paper  develops  the  mathe¬ 
matics  needed  to  describe  the  dynamic  interaction  of  the  frazil  ice  crystal  distribution  and  the  heat 
balance  of  the  turbulent  water. 

The  Literature  Review  section  covers  frazil  ice  in  natural  water  bodies,  and  industrial  crystalliza¬ 
tion. 

The  Basic  Lqitations  section  describes  the  dynamic  crystal  number  continuity  equation  in  parti¬ 
cle  phase  space  and  the  basic  equation  determining  the  heat  balance  of  the  water.  Also  defined  are 
the  moment  equations  defining  the  total  crystal  number,  area  and  volume. 


In  l  ho  lee  Crystal  (ire  wth  Runs  sod  ion.  liio  tii  .'\v  lit  late  along  llio  j-jxis  ( llio  m.ijoi  gi  ow  t  h  .imsI 
is  lound  to  bo  oontrollod  by  heal  liaiislei.  I  ho  ho.it  tiansloi  oooltloiout  toi  disk-shaped  panicles  is 
dotormmod  as  a  function  of  ciy stal  si/o  and  level  ot  iluid  luibuloiioo.  I  ho  magnitude'  >!  llio  moi Hal 
forces  and  buoyancy  forces  on  t ho  loo  oi>  sta Is  .no  dotoi iniuod  and  ilioii  inlluoiioo  oil  t ho  boat  iians- 
I or  assossod.  It  is  oonoludod  that  inoitial  toioos  can  bo  neglected  bin  that  giavitv  Imoo'  i u.i\  booomo 
significant  toi  large  crystals  in  weak  turbulence. 

Tito  \ ucliv n< ‘it  sootion  disoussos  both  initial  niioloation  and  sooond.uv  nuoloatioii  ot  Ii.i/iI  ll 
is  oonoludod  that  spontaneous  niioloation  ot  100  oan  bo  disoountod.  and  that  sood  oiystals  .no  noodod 
to  begin  frazil  growth  The  thooiotio.il  ratos  of  secondaiy  niioloation  aio  paitialk  mixleled.  and, 
t h o >  aro  found  to  dopond  on  tho  lovol  of  supercooling.  tho  turbulonoo  lo\ol  and  tho  ctyslal  si/o 
distribution. 

Tho  Irazil  lee  Dvnaimes  sootion  presents  a  simple  but  practical  anak  tioal  solution  ot  the  basic 
equations  tor  a  steady  -state  ory  stalli/or.  and  tho  results  are  presented  lot  tho  grow  tli  ot  Irazil  in 
such  a  orystalli/or.  At.  overview  ot  tho  basic  equations  :s  presented  emphasizing  the  intoiniatii 
feedback  inherent  in  tho  equations.  Tho  next  stops  in  future  research  are  outlined. 

LITERATURE  REVIEW 
Introduction 

Either  directly  or  indirectly,  tra/il  ice  has  drawn  the  attention  of  numerous  engineers  and  scien¬ 
tists,  and  the  literature  comes  from  many  fields.  Tho  first  part  of  this  section  is  a  brief  chronologi¬ 
cal  review  of  the  literature  produced  by  civil  and  mechanical  engineers,  who  are  concerned  with 
frazil  ice  in  northern  rivers  and  lakes,  the  meteorological  and  hy  draulic  conditions  tinder  which  it 
formed,  and  its  tremendous  negative  impact  on  man's  activity .  The  second  part  of  this  section  re¬ 
views  knowledge  gained  about  tra/il  trout  other  approaches,  primarily  from  interest  in  the  desalina¬ 
tion  of  sea  water  by  freezing.  In  this  context,  frazil  is  viewed  as  an  industrial  product  that  is  made 
in  bulk  crystallizers.  This  will  not  be  a  chronological  review,  but  rather  an  introduction  to  the  ex¬ 
tensive  literature  on  industrial  crystallization. 

Frazil  ice  in  natural  water  bodies 

Interesting  phenomena  such  as  ice  forming  on  fishing  lines  and  baskets  under  water,  lost  anchors 
being  lifted  by  masses  of  ice.  and  ice  displacing  buoys  by  moving  their  moorings  brought  early  atten¬ 
tion  to  frazil  ice.  But  traz.il  was  not  studied  seriously  until  people  began  to  use  rivers  and  lakes  un¬ 
interrupted  throughout  the  winter.  Hydropower,  winter  navigation  and  demands  for  fresh  watei 
by  factories  and  cities  provided  economic  incentives  for  the  study  of  frazil  ice. 

Research  on  tra/il  ice  in  natural  water  bodies  was  produced  in  response  to  several  major  concerns. 
The  first  was  predicting  where  frazil  ice  would  form.  The  ability  to  do  this  is  vitally  important  to 
those  who  must  keep  open  the  water  intakes  of  hydroelectric  plants,  water  supplies,  etc.  A  second 
was  predicting  the  total  volume  ol  tra/il  ice  that  a  given  reach  of  river  could  produce  because  large 
volumes  of  frazil  often  cause  dam-like  ice  jams  that  can  cause  catastrophic  Hooding.  A  third  major 
concern  was  coping  with  the  problems  caused  by  I  razil.  primarily  finding  ways  of  keeping  intakes 
operating,  keeping  canals  and  navigation  channels  open,  and  eliminating  ice  jams. 

The  meteorological  and  hydraulic  conditions  under  which  Irazil  can  form  were  determined 
through  observation  and  hard  experience.  Murphy  (  UK)1))  provided  a  good  description  of  these 
conditions.  Two  of  his  points,  that  Irazil  forms  in  areas  where  an  ice  cover  does  not  exist  and  that 
it  is  closely  associated  with  turbulent  waters,  are  universally  recognized. 

In  lee  I. 'ngincerin e.  Barites  ( l‘)2b)  summarized  the  knowledge  of  ice  formation  i.p  to  that  time, 
and  promulgated  certain  ideas  about  the  formation  of  frazil.  Many  ot  Barne-  ul,  were  not  cor¬ 
rect.  Unfortunately,  these  misconceptions  displayed  a  “remarkable  persistence"  (Varslens  l'tbb) 
and  many  years  passed  before  they  were  totally  rectified.  Barnes  did.  however,  firmly  establish  that 
there  was  always  supercooling  during  the  formation  of  Irazil. 


After  tin*  publication  of  Av  /  /»je///<vri»ri’  itt  I'-*"’  ..>q>ean  ifiVcMiga^is  ;--T.  ti.  .  .1  ! 

Altberg  in  Russia  and  Devik  in  \  jv  Ice  tosoaich  in  Kiismj  had  been  >:i:>  A.ac  :  ■  :  :  1. 

iee  completelv  *•'  ^i\eu  the  water  supply  of  the  eitv  01  Si .  I’eteishui e  in  1 1  *  1 4  ii.  ;li  I )  .  1  k  r  : 
Altberg  emphasized  quail  illative  accounting  ol  heal  lianstei  and  advanced  moie  ^  a  :  \  1  idea  •  ' ! 

mechanisms  of  ice  nucleauon  and  the  lole  ol  lonu-w ave  heal  transfei .  Both  believed  ' h a  1  '  -lei.a 
particles  suspended  in  the  water  piovuie  nuclei  tin  the  ice.  Both  emphasized  the  need  '  s;i;v: 
cooled  water  and  denionstialed  that  watei  bodies  pioducing  liazil  ice  weie  indeed  siip.-i. o. .!■  . i 
Dev  1  Is  ( Tarsi ei is  I ‘>00 )  is  ci edited  with  111 1 inducing  the  teiius  “active"  and  "passive"  1  o  disiu  .  .. ■  h 
frazil  ice  in  supercooled  water  (active)  liom  liazil  ice  in  watei  at  the  lieezine  point  (passive)  .i.v.. 
trazil  has  a  greater  tendency  to  stick  to  submerged  objects  than  passive  liazil.  \ltbeig  (  I'1  do  1  1:  le 
that  ice  crystals  introduced  into  supercooled  turbulent  watei  produced  a  great  uuiuhei  or  avldi’i  'ii. 
crystals,  and  suggested  that  portions  of  the  crystals  may  bieak  oil  and  serve  as  the  start  tot  a  new 
growth  ol  crystals.  This  idea,  which  has  been  expanded  into  the  tlieorv  ol  sccoiidatv  unclear i m. 
will  be  discussed  in  a  following  section . 

At  the  l‘>59  Montreal  Congress  of  the  International  Association  ol  llvdtaulic  Reseaich  ( lAIIKl. 
the  Committee  on  Ice  Problems  was  created.  Therefore,  the  ll>5(J's  can  serve  as  a  convenient.  1! 
arbitrary,  division  between  the  early  and  recent  liter  a  t  lire.  Seminars  on  ice  piohlems  weie  held  ai 
the  IAHR  Congresses  through  the  19(>0's.  In  ll>70  the  tirst  International  Svmposium  on  Ice  I’tob- 
lents  was  held  in  Iceland,  and  since  then  thev  have  been  held  regulatlv  ( J  1  l)"2 .  I 5 .  I 
1981  ),  Many  excellent  papers  on  frazil  are  contained  in  the  proceedings  of  these  symposia. 

Detailed  observations  of  frazil  formation  in  small  streams  were  published  by  Schaelet  (  1 1  < 5 < )  I 
and  Cilfillian  et  al.  ( 1972).  Wigle  ( 1970)  and  Arden  and  Wigle  ( 1972)  published  detailed  observa¬ 
tions  of  frazil  formation  in  the  Niagara  River.  Schaefer  observed  that  trazil  civ  stals  appealed  to 
form  initially  on  the  water's  surface  and  were  then  submerged  by  turbulence.  He  found  that  the 
frazil  crystals  retained  their  disk  shape  up  to  a  diameter  ol  approximately  2-7  mm  but  beyond  this 
dimension  dendrites  grew  out  from  the  flat  disks.  He  noted  that  these  dendrites  wet  •  quite  fragile 
and  easily  detached.  He  described  sponge-like  deposits  of  frazil  that  gather  011  the  upstream  side 
of  Tocks  and  other  underwater  objects,  but  he  drew  a  distinction  between  this  and  "true"  anchor 
ice  true  anchor  ice  being  sheet-like  crystals  of  ice  that  grew  out  from  submeiged  objects  to  which 
they  were  securely  fastened.  Schaefer  felt  that  true  anchor  ice  was  quite  rare  in  nature  and  that 
almost  all  reported  anchor  ice  was  actually  deposited  frazil. 

Gilfilhan  et  al.  ( 1972)  closely  observed  the  freezeup  of  a  small  subarctic  Alaskan  stream.  The 
border  along  the  edge  of  the  stream  grew  inwards  to  cover  it  5  days  after  freezeup  began.  They 
measured  an  increase  in  the  electrical  conductivity  of  the  stream  water  during  the  periods  of  Irazil 
production  and  they  attributed  this  to  the  rejection  of  impurities  by  the  ice  during  growth  of  the 
ice  crystals.  They  found  substantial  agreement  between  the  mass  of  Irazil  determined  by  estimat¬ 
ing  the  heat  loss  front  the  stream  and  the  mass  of  Trazil  determined  from  the  increase  in  electrical 
conductivity.  The  initial  ice  crystals  at  the  start  of  a  period  of  supercoolii  g  weie  discoids.  six- 
pointed  stars  or  hexagonal.  Trazil  crystals  larger  than  1  mm  in  diameter  had  scalloped  edges  or 
were  irregularly  shaped. 

Arden  and  Wigle  ( 1972)  observed  frazil  formation  on  the  upper  Niagara  Rivet .  a  laiee.  fast  (low  ¬ 
ing  river  that  remains  open  all  winter.  They  determined  a  “typical"  sequence  ot  ice  toimaiion  iluit 
took  place  each  night  when  meteorological  conditions  allowed  the  heal  content  ol  the  nvei  to  dis¬ 
sipate.  During  such  a  typical  sequence,  the  river  cooled  from  the  surface  down.  They  saw  eivstals 
on  the  surface  even  when  the  temperature  of  the  main  body  of  the  rivet  was  +0.07  (  .  As  1  ho  main 
body  ot  the  river  cooled,  the  crystals  on  the  surface  became  larger  and  were  submerged  to  greatei 
depths  by  the  turbulence.  Some  of  the  crystals  on  the  surface  weie  disks,  but  most  were  nregiilat . 
There  was  a  "transition"  period,  lasting  about  I  hour,  during  which  the  top  2-7  m  ot  water  was 
supercooled  while  the  bottom  5-7  m  retained  a  temperature  slightly  above  freezing.  I  razil  appeals' 
in  the  supercooled  layer  but  not  in  the  lower  layer;  however,  when  the  entire  1  ivei  became  supei 
cooled,  frazil  existed  throughout  the  depth  of  the  river.  During  this  period  disk-shaped  eivstals 
were  the  dominant  form. 


Figure  1.  Typical  supercooling  sequence  (after  Michel  I  ‘>6  ~ ). 


The  frazil  in  the  supercooled  water  immediately  stuck  tenaciously  to  any  submerged  object.  How¬ 
ever,  frazil  did  not  build  up  on  the  river  bottom  until  the  supercooled  water  had  existed  for  a  few 
hours.  During  a  typical  night  this  anchor  ice  buildup  reduced  the  flow  of  the  river  by  as  much  as 
25 7c. 

Beginning  in  the  1%0's  a  series  of  experiments  were  conducted  in  flumes  by  Michel  ( 1963)  and 
Carstens  (1966)  that  provided  very  consistent  qualitative  descriptions  of  frazil  formation.  Michel 
(1963)  determined  a  typical  sequence  of  water  temperatures  and  frazil  formation  (Fig.  I ).  At  the 
start  of  such  a  sequence  the  water  in  the  tlunte  cooled  at  a  constant  rate.  A  time  history  of  the 
water  temperature  formed  a  straight  line  with  constant  negative  slope  during  this  initial  period. 

When  the  water  temperature  reached  0°C,  it  passed  through  this  temperature  and  continued  cool¬ 
ing  at  the  same  rate  until  it  reached  a  few  hundredths  of  a  degree  below  the  freezing  point,  after 
which  the  rate  of  temperature  drop  began  to  decrease  rapidly  until  it  became  zero.  At  this  time 
tiny  ice  particles  formed  that  were  too  small  to  be  seen  by  the  naked  eye,  but  could  be  detected  if 
viewed  with  a  strong  light. 

After  the  water  reached  its  minimum  temperature,  it  returned  to  0°C.  Initially,  the  water  tem¬ 
perature  would  rise  quickly  and  then  more  and  more  slowly,  approaching  0°C  asymptotically.  The 
water  temperature  measured  by  Michel  during  the  course  of  an  experiment  reflected  a  balance  be¬ 
tween  the  heat  loss  at  the  water’s  surface  and  the  latent  heat  released  by  the  growing  frazil  crystals. 
When  the  rate  of  temperature  change  became  zero,  the  latent  heat  released  by  the  growing  crystals 
exactly  balanced  the  heat  loss  at  the  water  surface. 

Carstens  (1966)  conducted  a  series  of  experiments  in  a  test  flume  where  he  investigated  the  ex¬ 
tent  to  which  water  could  be  supercooled.  He  varied  the  rate  of  heat  loss  through  the  water  surface 
and  he  varied  the  turbulence  of  the  flow.  Carsten’s  experiments  produced  temperature-time  his¬ 
tory  curves  very  similar  to  those  of  Michel.  However,  the  temperature  of  the  water  did  not  always 
return  to  the  freezing  point',  it  would  often  reach  a  constant  value  of  supercooling,  which  he  termed 
the  residual  supercooling.  When  the  rate  of  heat  loss  from  the  water  surface  was  increased,  the  max¬ 
imum  supercooling  of  the  water  increased,  tire  residual  supercooling  increased  and  the  rise  in  tem¬ 
perature  from  the  maximum  supercooling  to  the  residual  supercooling  required  less  time.  Carstens 
also  felt  that  the  number  of  ice  crystals  increased.  He  saw  large  qualitative  differences  in  the  tem¬ 
perature-time  history  curves  when  the  turbulence  of  the  flow  was  varied.  “Strong"  turbulence 
produced  relatively  small  maximum  supercooling  levels,  and  the  temperature  would  return  rela¬ 
tively  quickly  to  the  freezing  point  with  little  or  no  residua!  supercooling,  while  “weak"  turbulence 
produced  relatively  large  maximum  supercooling  levels  that  would  persist  for  relatively  long  periods. 

Carstens  concluded  that  the  level  of  supercooling  reached  in  a  particular  experiment  reflected 
the  balance  of  the  heat  loss  at  the  surface  and  the  heat  transferred  from  the  growing  Ira/il  crystals. 
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which  was  controlled  by  turbulence  ot  the  flow .  lie  suggested  that  the  heat  tianstei  from  the  par¬ 
ticles  to  the  water  could  be  related  bs  an  expression  such  as  Nit  =  ('Re'"  where  Ntt  is  the  Nusselt 
number  (nondtmensional  heat  transfer).  C  and  in  are  physical  constants  in  Ins  nomenclatuie.  and 
Re  is  the  Reynolds  number  based  on  the  characteristic  velocity  of  the  (low  and  the  chataclctistic 
linear  dimension  of  the  I'ra/il  crystal. 

The  final  experiments  inspired  by  observations  of  Ira/tl  in  natural  water  bodies  to  be  described 
in  this  section  are  those  of  Muller  ( Il>78).  Muller  conuuncd  his  experiments  in  a  refrigerated  turbu¬ 
lence  jar  in  which  turbulence  was  generated  by  an  oscillating  grid.  I  le  investigated  the  rate  of  heat 
transfer  from  the  frazil  crystals  to  the  supercooled  water  and  the  rate  at  which  the  number  id'  frazil 
crystals  increased.  He  observed  no  frazil  at  supercooling  levels  up  to  I  C  unless  the  water  was  seed¬ 
ed  with  ice  crystals. 

The  number  of  ice  crystals  increased  during  the  course  of  an  experiment.  Muller  concluded  that 
some  type  of  multiplication  process  was  at  work,  but  he  could  gain  no  insight  into  the  mechanisms 
of  this  process.  However,  the  multiplication  process  “intensified”  with  an  increase  in  supercooling 
or  an  increase  in  turbulence.  He  found  that  the  apparent  heat  transfer  from  the  frazil  crystals  could 
be  normalized  with  the  supercooling  level  and  the  size  of  the  crystals  and  that  it  was  a  constant  for 
all  the  experiments  with  an  average  representative  value  of  the  Nusselt  number  of  10. 

Many  studies  have  attempted  to  predict  the  total  volume  of  frazil  that  can  be  produced  in  natural 
water  bodies  by  estimating  the  heat  budget  of  the  water  body.  Frevsteinsson  ( 1070)  and  Carstens 
( 1070)  presented  generalized  approaches.  Basically,  these  involved  determining  the  heat  budget  by 
accounting  for  the  total  heat  transfer  at  the  water  surface  by  evaporation,  convection  and  radiation, 
and  the  heat  transfer  from  other  sources,  such  as  through  the  bed  and  viscous  generation  of  heat. 
These  researchers  found  that  the  most  heat  was  transferred  through  the  water  surface,  and  that  the 
other  sources  could  generally  be  ignored.  They  assumed  that  once  the  water  temperature  reached 
the  freezing  point,  any  additional  heat  loss  resulted  in  the  production  of  frazil.  A  general  relation¬ 
ship.  M  =  QA/L,  was  used  to  calculate  A/,  the  mass  of  frazil,  where  /,  is  the  latent  heat  of  fusion, 

Q  the  total  heat  loss  per  unit  area,  and  A  the  surface  area  of  the  water  body. 

Estimates  of  the  amount  of  frazil  generated  in  the  Niagara  River  were  made  by  Ferguson  and 
Cork  ( 1972) and  List  and  Barrie  (1972).  Shen  ( 1980)  and  Slten  and  Ruggles  (1982)  made  very  de¬ 
tailed  estimates  of  the  volume  of  frazil  generated  in  the  open  sections  of  the  St.  Lawrence  River. 

In  all  these  studies  heat  losses  at  the  water  surface  dominated  the  heat  budget  of  the  water.  Shen 
found  that  the  bed  heat  flux  could  be  an  important  component, especially  for  colder  winters  when 
the  open  water  area  was  small. 

In  this  section  a  brief  review  of  the  civil  engineering  literature  relevant  to  frazil  has  been  presented. 
The  important  question  of  the  initial  nucleation  of  frazil,  which  has  been  long  debated,  will  be  dis¬ 
cussed  later  in  a  separate  section.  The  specific  problem  of  ice  control  at  water  intakes  was  not  men¬ 
tioned,  although  there  is  literature  on  the  subject  (Murphy  1 909,  Granbois  1953,  Logan  1974.  lan- 
tillo  1981 ).  In  addition,  several  very  good  general  reviews  of  frazil  have  been  published  ('Williams 
1959,  Michel  1971 .  Osterkamp  1978,  Ashton  1978,  Martin  1981). 

Frazil  ice  in  industrial  crystallizers 

Other  than  in  natural  water  bodies,  frazil  has  been  extensively  studied  only  in  the  industrial  crys¬ 
tallization  process  of  desalination  by  freezing.  Industrial  crystallization  is  the  production  of  cry  stals 
from  supersaturated  or  supercooled  solutions  in  agitated  vessels  of  varying  complexity  that  are 
called  crystallizers.  Generally,  the  average  crystal  size  and  crystal  size  distribution  are  very  impor¬ 
tant  in  determining  the  economic  returns  of  a  crystallizer  product.  Therefore,  investigators  have  ex¬ 
pended  much  time  and  effort  in  developing  both  theoretical  and  empirical  means  of  predicting  the 
crystal  size  distribution  that  will  result  front  a  given  crystallizer.  A  good  introduction  to  this  is 
Theory  <>J  Paniculate  Processes  (  Randolph  and  Larson  1971).  The  rigorous  determination  of  the 
crystal  distribution,  together  with  the  associated  heat  balance  and  the  appropriate  boundary  con¬ 
ditions,  provides  a  unified  predictive  and  descriptive  theory  for  desalination  by  freezing. 
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a.  Butane  freezing  process,  simplified  flow  diagram. 
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b.  Vacuum  flash  freezing  process,  simplified  flow  diagram. 
Figure  2.  Desalination  by  freezing. 


Desalination  by  freezing 

Desalination  by  freezing  has  been  extensively  studied  (Margolis  1969,  Kane  1971 ,  Hvans  1973. 
Simpson  et  a 1.  1973,  Woltz  1975,  and  Smith  and  Sarofirn  1979).  Desalination  by  freezing  is  pre¬ 
dicated  on  the  low  solubility  of  salt  in  ice.  The  ice  is  formed  in  an  agitated  crystallizer  from  which 
heat  is  removed;  the  supercooling  levels  achieved  in  the  crystallizers  are  small,  typically  less  than 
1°C.  The  dominant  shape  of  ice  formed  is  a  disk.  The  resulting  ice-brine  slurry  is  pumped  to  a 
washer  where  the  brine  is  drawn  off.  The  ice  is  then  melted  by  the  condensing  refrigerant. 

There  are  several  methods  for  removing  iieat  front  the  crystallizer.  The  secondary  refrigerant 
freezing  process  removes  heat  by  vaporizing  an  immiscible  refrigerant,  such  as  butane.  The  vacuum 
freezing-vapor  compression  process  vaporizes  a  portion  of  the  brine  at  low  pressure.  Both  processes 
are  shown  in  Figure  2. 

In  desalination  by  freezing,  the  quality  of  the  crystallizer  product  and  the  production  rate  largely 
determine  the  economic  returns  of  the  process.  The  product  quality  is  determined  by  the  average 
crystal  size  and  the  crystal  size  distribution,  and  it  is  important  in  determining  the  pressure  drop 
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across  lire  ice  bed  formed  in  the  washer  because  the  pressure  drop  is  inversely  piopoitumal  to  the 
square  of  the  crystal  diameter.  Also,  the  melting  rate  of  the  crystals  increases  as  the  crystal  M/e 
increases.  Therefore,  the  quality  is  improved  when  the  average  si/e  is  increased  and  when  the  num¬ 
ber  of  small  particles  is  decreased.  The  production  rate  is  equal  to  the  crystalli/er  volume  times  the 
slurry  density  divided  bv  the  average  residence  time  ol  the  ice  in  the  crystalli/er.  Thus  the  produc¬ 
tion  rate  can  be  improved  if  the  slurry  density  is  increased,  it  the  average  crystal  si/e  is  increased  or 
if  the  residence  time  is  decreased. 

Generally,  the  type  of  crystalli/er  proposed  for  this  process  would  operate  at  steady  state  and  is 
called  a  Mixed  Suspension,  Mixed  Product  Removal  (MSMPR)  crystalli/er.  Operation  at  steady  state 
requires  that  the  number  of  crystals  and  the  crystal  si/e  distribution  do  not  change  with  time.  New 
crystals  must  be  nucleated  in  the  crystalli/er  as  the  crystallizer  is  operated  with  clear  (unseeded)  feed 
streams.  At  steady  state,  one  crystal  must  be  nucleated  for  each  crystal  in  the  crystalli/er  in  an  aver¬ 
age  of  one  residence  time.  Therefore,  the  residence  time  is  equal  to  the  number  of  crystals  (.V)  di¬ 
vided  by  the  nueleation  rate  present  in  the  crystalli/er  (,V).  Tire  average  crystal  size  can  be  deter¬ 
mined  by  multiplying  a  representative  growth  rate  (O')  of  the  crystals  by  the  average  residence  time. 
The  average  crystal  size  ( r )  is  then 

7  =  GN/N. 

The  average  crystal  size  is  a  function  of  the  nueleation  and  growth  rates.  The  nueleation  rate  and 
the  growth  rate  are  therefore  very  important  in  determining  the  economics  of  desalination  by  freezing. 

Nueleation  rate 

The  type  of  nueleation  in  ice  crystallizers  is  secondary  nueleation,  and  secondary  nueleation  takes 
place,  irrespective  of  its  mechanism,  only  because  of  the  presence  of  crystals  of  the  material  being 
crystallized  (Botsaris  1976).  For  secondary  nueleation  to  occur,  the  existing  stable  crystals  must  in 
some  way  generate  new  crystals.  Inquiry  into  this  process  can  be  conveniently  divided  into  two 
questions:  what  is  the  source  of  the  potential  nuclei  of  the  new  crystals  and  how  are  these  nuclei 
removed  from  the  source  and  displaced  into  bulk  solution  to  initiate  new  crystals? 

A  number  of  sources  of  the  ice  nuclei  have  been  suggested,  including  dendrites,  microscopic  sur¬ 
face  irregularities  and  large-scale  breakage  of  the  parent  crystal.  Dendrites  and  surface  protuberances 
are  often  reported  on  frazil  crystals  and,  if  detached  from  the  parent  crystal,  could  easily  initiate 
new  crystals.  However,  there  is  ample  evidence  that  secondary  nueleation  can  happen  without 
visible  dendritic  growth  (Evans  1973,Woltz  1975).  Clontz  and  McCabe  (1971 )  were  able  to  generate 
crystal  nuclei  by  causing  crystals  to  collide  with  rods  of  various  materials.  The  contacts  produced  no 
visible  defects  in  the  parent  crystals,  but  they  assumed  that  during  the  collisions  microscopic  surface 
irregularities  were  sheared  off.  These  microscopic  particles  are  considered  the  most  likely  source  of 
potential  ice  crystal  nuclei  (Evans  et  a!.  1974a,  b).  Large-scale  breakage  of  ice  crystals  is  rarely  ob¬ 
served. 

The  production  of  the  crystal  nuclei  from  the  sources  described  above  has  been  attributed  to  fluid 
shear  and  collisions  of  all  types.  Clontz  and  McCabe  ( 1971 )  found  no  evidence  that  fluid  shear  could 
produce  new  crystals.  Secondary  nueleation  induced  by  collisions  among  the  parent  crystals  has  been 
investigated  by  Lai  et  al.  (1969),  Garabedian  and  Strickland-Constable  (1972).  Ottens  el  al.  ( 1972), 
Denk  and  Botsaris  (1972)  and  Desai  et  al.  (1974).  Secondary  nueleation  of  ice  crystals  by  collision 
with  other  surfaces  has  been  investigated  by  Garabedian  and  Strickland-Constable  ( 1 974),  Evans 
et  al.  ( 1974a.  b),  and  Woltz.  (1975).  The  number  of  crystals  produced  per  collision  was  found  to 
increase  with  the  collision  energy  and  the  supercooling  of  the  fluid . 

Garabedian  and  Strickland-Constable  ( 1974)  proposed  the  survival  theory  to  explain  these  re¬ 
sults.  According  to  the  survival  theory,  the  number  of  nuclei  produced  per  collision  is  proportional 
to  the  collision  energy  and  independent  of  the  supercooling  level.  The  nuclei  produced  are  distribut¬ 
ed  over  a  range  of  sizes,  and  only  those  larger  than  a  critical  radius,  determined  by  the  supercooling 
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of  the  solution,  survive  and  grow.  Recent  observations  ((iarside  and  I  .at  son  I97N)  have  produced 
doubts  about  some  of  the  details  of  the  survival  theory;  howevei .  the  models  ol  secondaiv  nuclea- 
tion  based  on  its  premises  are  very  general  and  can  incorporate  these  results. 

The  rate  of  production  of  new  crystals  by  secondary  nueleation  can  be  limited  by  eitliei  the  pro 
duction  rate  of  the  potential  nuclei  on  the  crystal  surface  or  the  rate  of  their  removal  tiom  i he  cry  s¬ 
tal  surface.  Being  limited  by  the  removal  rate  implies  that  the  time  between  collisions  is  long  enough 
to  allow  the  crystal  surface  to  heal  from  previous  collisions,  so  the  nueleation  late  is  therefore  inde¬ 
pendent  of  the  system  history.  The  secondary  nueleation  of  ice  is  limited  by  the  removal  talc, 
bvans  et  al.  ( 1974b)  proposed  that  the  overall  nueleation  rate  (A'r)  for  ice  with  more  than  one 
mechanism  of  removal  is  the  linear  sum  of  the  actual  nueleation  rale  altiibutable  to  each  mechan¬ 
ism  of  removal.  Thus 

,VT  =  ,V1+A’,  +  /V,3+...;Vl 

where  AA  is  the  nueleation  rate  attributable  to  the  / th  mechanism. 

For  example,  in  a  crystallizer,  the  overall  nueleation  rate  is  the  linear  sum  of  the  nueleation  rate 
caused  by  collisions  of  the  ice  crystals  with  the  agitator  (propelior  or  turbine),  the  cry  stalli/er  walls 
and  other  crystals.  Evans  et  al.  ( 1974b)  demonstrated  that  coating  the  agitator  and  crystalli/er 
walls  with  a  material  that  is  softer  than  ice  could  reduce  the  nueleation  rate.  The  coating  reduced 
the  contact  energy,  thereby  reducing  the  number  of  nuclei  produced  per  collision. 

Crystal  growth  rates 

The  growth  of  crystals  in  agitated  turbulent  solutions  has  long  been  a  topic  of  research  by  chemi¬ 
cal  engineers.  An  extensive  review  is  provided  by  Wadia  ( ll)74).  The  growth  mechanisms  for  ice 
crystals  suspended  in  supercooled  water  are  the  incorporation  of  water  molecules  into  the  ice  sur¬ 
face  (the  crystalline  kinetics)  and  the  transport  of  the  latent  heat  of  fusion  to  the  bulk  fluid.  Thus, 
the  overall  growth  rate  of  an  ice  crystal,  in  fresh  or  salt  water,  is  the  result  of  the  interaction  of  the 
crystalline  kinetics  and  the  transport  processes  of  the  surrounding  turbulent  fluid.  It  has  been 
found  that  the  growth  rates  of  the  major  axis  of  disk  crystals  are  largely  controlled  by  transport 
processes.  The  transport  processes  can  be  described  if  the  geometry  of  the  crystal  is  known  and  if 
the  ambient  velocity  distribution  of  the  fluid  in  the  vicinity  of  the  crystal  can  be  adequately  de¬ 
scribed. 

Unfortunately,  the  fluid  in  a  highly  agitated  crystallizer  exhibits  a  turbulent  How  Held  with  many 
seemingly  random  aspects.  However,  successful  transport  correlations  predicated  on  Kolmogomv's 
theory  of  isotropic  turbulence  (Hinze  1 959)  have  been  developed.  In  a  turbulent  crystallizer,  the 
impeller  adds  a  known  amount  of  energy  and  creates  isotropic  turbulence  around  the  suspended 
crystals.  The  ambient  velocity  distributions  are  then  determined  by  the  energy  put  in  by  the  im¬ 
peller,  the  viscosity  of  the  fluid  and  the  crystal  diameter.  The  transfer  processes  are  calculated  by 
determining  the  heat  transfer  given  the  ambient  velocity  distribution.  The  specific  case  of  heat 
transfer  from  suspended  frazil  crystals  will  be  described  in  more  detail  in  the  Ice  Crystal  Growth 
Rates  section. 

Summary 

In  summary,  desalination  by  freezing  is  a  process  of  industrial  crystallization  in  which  frazil  ice 
is  purposefully  created  to  make  fresh  water  from  salt  water.  To  provide  the  maximum  economic 
return,  the  crystal  size  distribution  has  to  be  optimized  as  described.  The  two  major  influences  on 
the  crystal  size  distribution  are  the  growth  rate  and  the  rate  of  secondary  nueleation  ot  the  ice 
crystals,  and  as  a  result,  the  basic  mechanisms  that  control  these  rates  have  been  extensively  Re¬ 
searched.  This  research  has  not  yet  been  applied  to  frazil  production  in  natural  water  bodies 


BASIC  EQUATIONS 


Introduction 

In  this  section  a  descriptive  and  predictive  crystal  distribution  theory,  following  Randolph  and 
Larson  (1971),  will  be  developed.  Rigorous  determination  of  the  population  balance  of  crystals 
will  result  in  a  crystal  number  continuity  equation,  the  solution  of  which  is  a  distribution  function 
that  contains  information  about  the  form  and  the  magnitude  of  the  crystal  distribution.  The  func¬ 
tion  describes  the  population  density  of  crystals  of  each  specific  linear  dimension  throughout  the 
volume.  The  physical  parameters  that  affect  the  crystal  size  distribution  appear  in  the  crystal  num¬ 
ber  continuity  equation.  An  additional  equation  that  describes  the  heat  balance  within  the  vol¬ 
ume  is  required  to  determine  the  crystal  size  distribution.  The  two  equations  are  linked  by  the 
growth  and  secondary  nucleation  rates  of  the  ice  crystals,  which  are  dependent  on  both  the  heat 
balance  and  the  crystal  size  distribution.  In  theory  the  equations  can  be  solved  if  the  various  re¬ 
quired  boundary  and  initial  conditions  are  known.  These  equations  will  serve  as  the  framework  for 
this  investigation  of  frazil  ice  in  natural  water  bodies.  In  addition,  the  crystal  moment  equations 
will  be  developed.  These  equations  are  very  convenient  for  determining  the  total  mass  and  surface 
area  of  the  crystals  as  a  function  of  size. 

Crystal  number  continuity  equation 

The  crystal  distribution  will  be  described  in  a  space  termed  the  crystal  phase  space  or  more  gen¬ 
erally  the  particle  phase  space.  Particle  phase  space  is  defined  by  the  least  number  of  independent 
coordinates  that  provides  a  complete  and  useful  description  of  the  properties  of  the  crystal  distribu¬ 
tion.  It  is  convenient,  if  somewhat  arbitrary,  to  divide  particle  phase  space  into  two  subregions  de¬ 
fined  by  external  coordinates  and  internal  coordinates.  The  external  coordinates  describe  the  spa¬ 
tial  distribution  of  the  crystals.  Internal  coordinates  refer  to  properties  attached  to  each  individual 
crystal,  which  quantitatively  measure  its  state,  and  are  independent  of  its  position. 

To  begin,  a  crystal  distribution  function  n(R,  r )  will  be  considered.  This  function  is  defined 
over  a  region  R  of  the  particle  phase  space  consisting  of  the  three  spatial  dimensions  (the  external 
coordinates)  plus  any  number  of  internal  property  coordinates.  In  all  further  cases  the  internal  co¬ 
ordinates  will  be  restricted  to  one,  which  will  correspond  to  a  major  linear  dimension,  r,  of  the  ice 
crystals.  The  function  n(R,  t)  is  defined  as  the  population  density  of  crystals  in  the  region  R.  At 
a  time  t  the  number  of  crystals  in  an  incremental  region  of  the  particle  phase  space  JR  is  given  by 

dN  =  ndR  ( 1 ) 

and  the  total  number  of  crystals  in  a  region  R  at  time  t  is 

N(R)  =  / ndR.  (2) 

k 

Individual  crystals  can  continuously  change  their  position  in  the  particle  phase  space.  If  these 
changes  are  regular,  that  is,  if  they  are  the  result  of  gradual  and  continuous  movement,  the  convec¬ 
tive  crystal  velocity  along  a  respective  particle  coordinate  can  be  defined.  Let  I'e  be  the  external 
convective  velocity  and  be  the  internal  convective  velocity.  I'e  is  the  velocity  of  the  crystal 
through  space  and  will  result  from  the  interactions  of  all  the  forces  acting  on  the  crystal.  1',,  as 
restrictively  defined  here,  simply  represents  the  rate  of  change  of  the  size  of  a  crystal.  The  veloc¬ 
ity  of  the  overall  vector  crystal  phase  space  is  defined  as 

nR.t)=VAR.ft+VAR,t).  (3) 
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li  may  also  be  necessary  to  deal  with  the  sudden  appearance  (birth)  or  disappearance  (death)  of 
crystals  at  a  point  in  the  particle  phase  space.  For  example,  the  breakage  of  a  crystal  would  cause 
the  sudden  disappearance  of  a  crystal  at  the  internal  coordinate  corresponding  to  its  si/e  and  the 
sudden  appearance  of  its  fragments  at  the  coordinates  corresponding  to  their  si/es.  The  net  appear¬ 
ance  in  an  incremental  region  dR  at  a  time  would  be 

(B-D)dR  (4) 

where  B{R.  t)  and  D(R.  t )  represent  birth  and  death  functions  at  a  point  in  the  phase  space. 

The  population  balance  of  crystals  in  some  fixed  region  /f,  which  moves  convcctivcly  with  the 
particle  phase  space  velocity  F,  can  be  defined  as 

~  I  fu/R  =  I (B-D)ilR  .  (5) 

at  ■'  ■ 

R  K 

Expanding  the  first  term  using  Leibnitz’s  rule,  and  noting  that  the  region  R  was  arbitrary,  we  can 
state  the  population  balance  in  general  terms  as 

+V(Fe'0+V(Fj//)-S+Z)  =  0  .  (6) 

Tltis  is  the  number  continuity  equation  in  particle  phase  space.  This  equation  is  quite  general,  as 
each  term  is  defined  over  the  three  spatial  coordinates,  the  internal  coordinate  r  and  time.  As  only 
a  single  internal  coordinate  r  is  considered,  the  divergence  of  the  internal  convective  velocity  can 
be  represented  as 


v(fV0  =  iL  (G'/t) 

or 


(7) 


where  G(r.t)  is  the  convective  velocity  along  r  or  simply  the  growth  rate  of  the  ice  crystal  such  that 


G  = 


dr 
dt  ' 


(8) 


Equation  6  is  then 

~  +  ~(Gn)+D-B  +  V(Ven)  =  0. 


(9) 


This  is  the  number  continuity  equation  in  general  form.  Further  simplifications  can  be  made  by 
applying  physical  knowledge  of  frazil  and  the  conditions  under  which  it  evolves.  In  principle,  eq  9 
can  be  solved  for  n  (crystal  size  distribution)  if  G,  B  and  D  are  known  functions,  and  the  initial  and 
boundary  conditions  of//  are  known.  To  uniquely  define  these  parameters,  the  heat  balance  of  the 
volume  must  be  known.  The  basic  equation  governing  the  heat  balance  will  be  developed  later. 

Crystal  moment  equations 

It  is  often  of  interest  to  know  the  distribution  with  crystal  size  of  the  number  of  crystals,  the 
mass  of  crystals  or  surface  area  of  the  crystals.  These  can  be  uniquely  determined  if  the  distribu¬ 
tion  function  is  known.  Let  n(r,  t)  be  the  distribution  of  the  crystal  population  along  the  crystal 
size  axis  rat  a  time  t.  For  convenience  we  assume  that  n  is  not  dependent  on  any  other  coordinate 
and  t  is  fixed;  n(r)  is  the  population  density  of  crystals  of  size  r.  By  definition,  the  concentration 
of  crystals  of  size  r  to  r+dr  found  in  a  unit  volume  is 
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di\'  =  »(r)Jr. 


(10) 


The  units  ot'  n  are  (r-3 ),  (r  1 )  or  (r"4 ).  Tire  total  concentration  of  cry  stals  is 

;V  =  j  n(r)dr.  (II) 

o 

The  units  of  A'  are  crystals  per  unit  volume.  The  surface  area  of  an  individual  cry  sial  can  be  repre¬ 
sented  as 


a(r)  =  Kar2  (12) 

where  Ka  is  a  shape  factor  relating  area  to  size  squared.  For  geometrically  similar  crystals  the  shape 
factor  is  independent  of  size.  The  total  surface  area  per  unit  volume  of  suspension  is 

Ar  =  A'a  j  r 2  n(r)dr.  (Id) 

o 

The  weight  of  an  individual  crystal  can  be  represented  as 

m(r)=piKvr3  (14) 

where  ATV  is  a  volumetric  shape  factor  relating  the  volume  to  size  cubed  and  p  is  the  crystal  density. 
The  total  mass  of  crystals  per  unit  volume  of  suspension  is 

Mr  -  pjXv  j  r 3  n(r)dr  .  (15) 

o 

Heat  balance 

In  this  section  the  general  expression  for  the  heat  balance  of  the  frazil  ice-water  system  will  be 
developed.  This  expression  will  be  developed  strictly  for  frazd  crystals  suspended  in  fresh  water. 

Consider  a  differential  volume  in  which  pf  is  the  mass  concentration  of  water  (grams  of  water 
per  cubic  centimetre  of  mixture)  and  pt  is  the  mass  concentration  of  ice.  The  temperature  of  the 
water  is  7}  and  the  ice  Tke,  then 


C  A  (pf  Tf)  +  Cj  A  (p(  7'jce)  +  Cv(pf  Tf  Pf)  +  C,  V(Pj  Tke  fe )  = 


-  VA  vrf  +  pfu<h  +  Q*  +  /.fe  A  p,+v(Pj  p^j 


( lba) 


where  C=  heat  capacity  of  liquid  water 
Cj  =  heat  capacity  of  ice 
Ff  =  convective  velocity  of  the  fluid 
k  =  thermal  conductivity 
u  =  kinematic  viscosity 
d>  =  dissipation  function  (viscous  heating) 

Q*  =  net  heat  transfer  from  the  mixture 

Lte  =  latent  heat  of  fusion  of  the  ice  at  the  equilibrium  temperature  of  (he  mixture. 


Several  simplifications  can  be  made  to  cq  16a.  The  first  is  to  assume  that  7'f  '  7’K(.,  thai  is.  the 
temperature  of  the  ice  is  that  of  the  water.  While  this  cannot  be  strictly  true  (if  it  were  there  could 


be  no  heat  transter  between  the  ice  and  water),  this  assumption  introduces  onk  a  very  small  eiroi. 
Tile  error  is  small  as  ihe  temperatures  of  the  iee  and  water  are  very  elose.  iisiuih  within  0.1  ('.  AIm 
the  mass  of  ice  is  usually  very  small  compared  to  the  mass  of  water  (  *  I  ).  so  that  to  good  approxi¬ 
mation.  pf  can  be  set  equal  to  pf.  the  density  of  water,  and  p,  +  P,  ^  Pt  '  P,  .  A  second  assumption 
is  then  C\  P-JC  pf  «  I .  Additionally,  heat  conduction  can  he  neglected  and  heal  capacities  and  the 
latent  heat  can  be  considered  constant  because  of  the  small  variation  in  temperature.  Therefore, 
based  on  the  above  assumptions,  eq  I  Pa  can  be  rewritten 


dTf  -  / 

■F+V(r'r^'c^ 


d  -*•  O* 

—  p  +V(  r  p  )  +  — —  + 

a  t  K  eP'  cpf  c 


u‘l> 


(lob* 


where  L  is  the  latent  heat,  which  can  be  assumed  to  be  a  constant. 


Parameters  in  the  basic  equations 

Tire  two  basic  equations  are  the  crystal  number  continuity  (eq  9),  and  the  heat  balance  (eq  W>bl- 
The  various  parameters  that  appear  in  these  two  equations  will  now  be  discussed. 

1 .  G-The  growth  rate  of  the  major  linear  dimension  of  the  crystals  is  a  function  of  the  heat  trails 
fer  and  the  intrinsic  kinematics  of  the  ice  crystal.  In  the  Ice  Crystal  Growth  Rates  section  it  will  be 
shown  that  G  is  effectively  determined  by  the  heat  transfer  rate.  Thus,  in  general 

G  =  !liLJl  q 
P  -,  L 

where  h  is  the  heat  transfer  coefficient  and  is  a  function  of  the  crystal  size  r  and  the  level  of  turbu¬ 
lence  e:  0  is  the  supercooling  of  the  mixture. 

2.  D-The  death  function  can  be  set  to  'ero  for  all  signs  of  crystals.  This  is  equivalent  to  assum¬ 
ing  that  there  is  no  large-scale  breakage  of  the  crystals. 

3.  B-  The  birth  function  is  determined  by  the  rate  of  the  sudden  appearance  of  new  crystals. 

New  crystals  can  appear  as  a  result  of  spontaneous  nucleation.  secondary  nucleation  and  the  intro¬ 
duction  of  crystals.  The  Nucleation  section  will  show  that  spontaneous  nucleation  is  not  possible 
under  frazil-forming  conditions.  Therefore  B  will  be  determined  by  the  rate  at  which  new  crystals 
are  introduced  and  the  rate  of  secondary  nucleation. 

Let  iV-j-be  the  rate  of  secondary  nucleation.  /VT  is  a  function  of  the  crystal  distribution  n,  the 
turbulence  dissipation  rate  e,  the  supercooling  of  the  mixture  0.  and  perhaps  other  parameters.  Let 
Nt  be  the  rate  at  which  new  crystals  are  introduced.  We  assume  that  new  crystals  arc  created  and 
introduced  at  a  size  small  enough  that  the  radius  of  new  crystals  can  be  approximated  as  zero.  Thus 

fl  =  [yvT(0.  n,  e)+yV,J  5(r-0) 

where  6(r-0)  is  the  dirac  delta  function  [5(r  =  0)  =  1 , 6(r  ^  0)  =  0 1 . 

4.  Pj-The  mass  of  ice  per  unit  volume  of  mixture  can  conveniently  be  determined  using  the 
moment  equation 

P,  =  P,KV  I'  r'  n(r)dr  . 
o 

5.  Q*.  <1>-These  functions  are  determined  by  the  environment  of  the  water  body  of  interest, 
and,  in  particular,  the  meteorologic  and  hydraulic  conditions.  '1>,  the  dissipation  function,  will  gen¬ 
erally  be  quite  small  and  under  most  circumstances  can  be  set  equal  to  zero  with  small  error,  re¬ 
pressions  to  determine  the  va'ue  of  these  functions  will  not  be  developed  in  this  report. 

6.  I’  ,  lf  The  convective  velocity  of  the  ice  crystals  and  the  fluid  will  gcnerallv  be  very  similar. 


The  action  of  buoyancy  and  inertial  forces  on  the  ice  crystals  may  cause  the  ice  crystal  velncits  to 
differ  from  the  fluid  velocity  if  these  forces  become  large  compared  to  the  fluid  drag  force. 
Substituting  the  above  expressions  into  eq  9  and  16b  gives  us 

^  ^  ( h  n )  +  I ; 'n )  =  (TV,  +  A', )  5<r  -  0 ) 

dr  pL  dr  e  11 


dO 


L  Pi  I  a 


Q* 

Cpt 


where  n  =  n(x,y.z,r,t) 

0  =  0(x,y.z,t) 
h  =  hy,e) 
i;  =  lyix.y.z.r.t) 

A'r  =  .Vj(fl.e.M) 

If  =  If(.v,.v,r,r) 
Q*  =  Q*{x,y.z,t ). 


Writing  the  equations  in  this  form  emphasizes  the  dynamic  way  in  which  they  interact.  To  deter¬ 
mine  0  and  n  uniquely,  both  equations  must  be  solved  simultaneously,  and  the  boundary  condi¬ 
tions  and  initial  conditions  of  0  and  n  must  be  known.  Difficulties  arise  because  0  and  n  are  di¬ 
mensionally  incompatible.  In  the  Ice  Crystal  Growth  Rates  section  expressions  for  h  will  be  devel¬ 
oped  and  in  the  Sucleation  section  expressions  for  ,Vr  will  be  developed. 


ICE  CRYSTAL  GROWTH  RATES 
Introduction 

This  section  is  a  discussion  of  all  the  relevant  factors  that  determine  the  growth  rates  of  ice 
crystals  suspended  in  turbulent  flows.  These  factors  include  the  morphology,  the  intrinsic  kinetic 
growth  rates  and  the  heat  transfer  processes.  The  overall  goal  of  this  section  is  to  determine  the 
value  of  the  parameter  G  that  appears  in  the  number  continuity  equation. 

Morphology 

A  description  of  the  morphology  of  ice  is  not  simple.  The  various  shapes  of  ice  crystals  appear 
to  result  from  a  complex  interaction  of  the  imposed  heat  transfer  conditions  and  the  intrinsic  crys¬ 
tallography  of  ice.  We  know  from  observations  that  the  dominant  shape  of  ice  crystals  that  grow 
at  the  supercooling  levels  found  in  turbulent  water  bodies  is  a  flat  disk.  Virluallv  all  field  observa¬ 
tions  of  frazil  ice  note  that  the  crystals  are  disk  shaped.  It  has  been  reported  that  ice  crystals  in  the 
shape  of  six-pointed  stars,  hexagonal  plates  or  spheres,  and  small  pieces  of  dendritic  ice  all  evolve 
into  the  disk  shape  in  natural  water  bodies.  Researchers  have  studied  the  morphology  of  large 
numbers  of  frazil  ice  crystals  because  of  interest  in  desalination  by  freezing.  The  observations  of 
Margolis  ( 1969)  indicate  that  the  thickness  of  the  frazil  disk  is  0.68r  +  16.7'/?,  where  r  is  the  major 
radius  of  the  disk.  Smith  and  Sarofim  ( 1979)  say  that  t he  maximum  radius  is  approximately  0.8 
mm  for  disk  crystals  produced  in  turbulent  crystallizers. 

Disk-shaped  crystals  have  been  studied  in  the  laboratory  by  Kumai  arid  Itagaki  (1953),  Arakawa 
( 1954)  and  Williamson  and  Chalmers  ( 1966).  Arakawa "s  classic  drawing  of  the  growth  of  a  disk 
crystal  is  shown  in  Figure  3.  Arakawa  created  disk  crystals  by  first  growing  dendritic  ice  crystals 
on  the  bottom  of  a  container  containing  slightly  supercooled  water.  He  then  scratched  the  crystals 
with  the  tip  of  a  line  needle.  Spherical  ice  particles  with  a  diameter  ol  about  10  2  mm  formed 
and  rose  towards  the  surface.  Two  flat  spots  formed  on  the  spherical  particles’  surtace  as  they 


floated  up.  and  further  growth  was  lateral 
the  pat  tides  became  disks.  The  crystals  grew 
into  disks  with  diameters  of  0.5  to  5  nun  and 
with  diameter-to-thickness  ratios  of  5 : 1  to 
100: 1 .  Once  the  disks  grew  to  a  certain  st/e 
they  took  on  a  notched  look  with  many  notch- 
shaped  growths  at  the  perimeter  of  the  disk. 

In  a  series  of  experiments,  Williamson  and  Chalmers  ( 1060)  found  that  at  supercooling  levels  of 
0.2°C  or  less,  ice  grows  into  disk  crystals.  On  increasing  the  supercooling  level  to  0.4T,  they  found 
that  the  disk  morphology  became  unstable  and  small  protuberances  appeared  on  the  edge  of  the 
disks.  The  protuberances  did  not  appear  to  show  any  preferred  growth  directions.  At  supercooling 
levels  greater  than  approximately  0.6°C.  the  ice  assumed  the  characteristic  dendritic  appearance 
shown  by  many  snowflakes,  and  preferred  growth  directions  developed.  For  supercooling  greater 
than  -1°C,  preferred  growth  directions  dominated  the  morphology  and  produced  dendritic  crystals 
with  at  least  three  generations  of  branches. 

Researchers  think  that  the  disk  shape  taken  by  ice  at  low  supercooling  levels  is  the  result  of  the 
anisotropic  growth  kinetics  of  ice.  Ice  has  two  growth  axes  and  they  are  identifiable  by  their  opti¬ 
cal  properties.  The  ice  molecule  is  hexagonal.  The  hexagonal  axis  is  the  r  -axis.  and  the  three  axes 
normal  to  the  c-axis  are  the  a-axes.  The  a-axes  are  all  equivalent  because  the  crystal  is  symmetrical 
about  the  c-axis.  The  plane  that  contains  the  a-axes  is  the  basal  plane,  and  the  growth  rate  in  it  is 
tens  of  times  faster  than  that  parallel  to  the  c-axis. 

The  disk  shape  of  the  ice  crystal  is  subject  to  morphological  instability.  In  natural  water  bodies 
this  instability  is  observed  as  scalloped  edges,  dendritic  growths,  irregular  protuberances,  etc.,  on 
the  perimeter  of  the  crystals.  Experiments  have  shown  that  the  disk  shape  always  becomes  unstable 
if  the  supercooling  is  increased  enough  or  if  the  diameter  of  the  crystal  becomes  large  enough.  In¬ 
stability  is  perhaps  the  mechanism  that  limits  the  maximum  disk  radius  to  0.8  mm.  The  entire 
reason  for  the  instability  of  the  disk  shape  is  not  known,  although  it  has  been  the  subject  of  several 
investigations  (Williamson  and  Chalmers  1 966,  Fujioka  and  Sekerka  1 974).  Williamson  and  Chal¬ 
mers  concluded  that  the  instability  of  the  disk  shape  depends  on  heat  How  into  the  liquid,  not 
crystallography.  This  conclusion  seems  correct,  as  the  diameter  at  which  the  crystal  becomes  un¬ 
stable  varies  between  experiments. 

Intrinsic  kinetic  growth  rate 

The  mechanisms  that  determine  the  rate  at  which  an  ice  crystal  can  grow  are  transport  of  water 
molecules  to  the  crystal  surface  (for  ice  grown  in  the  pure  water  this  is  not  a  consideration),  their 
incorporation  into  the  crystal  surface  and  the  transport  of  latent  heat  away  from  the  surface.  The 
incorporation  of  molecules  is  controlled  by  the  crystallization  or  interface  kinetics  of  ice  and  the 
heat  transfer  reflects  the  particular  physical  situation  of  the  system  under  consideration.  The  inter¬ 
face  kinetics  of  ice  has  been  studied  both  theoretically  and  experimentally.  As  noted  in  the  previ¬ 
ous  section,  ice  has  two  principal  growth  directions:  within  the  basal  plane  (a-axis  growth)  and 
normal  to  the  basal  plane  (c-axis  growth).  The  interface  kinetics  of  each  growth  direction  appear 
to  be  different  and  each  will  be  addressed  separately. 

The  crystallization  kinetics  of  ice  can  be  viewed  as  part  of  the  larger  study  of  crystal  growth 
from  a  melt,  which  has  developed  an  extensive  literature.  In  the  interest  of  brevity,  the  interface 
kinetics  will  be  divided  into  three  simple  possibilities,  following  the  example  of  Fletcher  ( 1970). 

The  first  possibility  is  that  of  a  perfectly  smooth  crystal  face.  Because  of  the  lack  of  near  neigh¬ 
bors,  it  is  relatively  difficult  for  molecules  to  become  attached  to  such  a  surface.  It  is  necessary  lor 
a  stable  “island”  to  nucleate  on  the  surface  betore  growth  can  begin.  For  this  type  of  growth  it 
can  be  shown  that 
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Figure  . 1.  Frazil  disks  (after  Arakawa  I9>4). 


I IX) 


<7k  *  exp  (-a/0,) 


14 


where  G'k  is  the  kinetic  growth  velocity,  a  is  ;i  constant  and  u,  is  the  supercooling  ai  the  iniei  lace. 
This  type  of  growth  is  called  surface  nucleation. 

Real  crystal's  surfaces  are  generally  not  perfectly  smooth  and  may  contain  imperfections  of  van- 
ous  kinds.  A  dislocation  imperfection  allows  the  cry  stal  surface  to  advance  continuoush  .  I  01  this 
type  of  growth  it  can  be  shown  that 

(l‘M 

This  type  is  called  grow  th  by  dislocation  mechanisms. 

A  crystal  surface  may,  under  certain  circumstances,  be  rough  on  a  molecular  scale.  AI!  crystals 
show  a  general  tendency  to  become  rough  at  large  supercooling  levels  and  it  is  relatively  easy  tor 
molecules  to  become  attached  to  a  rough  surface.  The  velocity  of  this  type  of  grow  th  can  he  shown 
to  be 
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This  type  ol  growth  is  called  continuous  growth.  The  continuous  growth  mechanism  may  be  al- 
tected  by  the  curvature  of  the  interface,  especially  if  the  radius  of  curvature  is  very  small. 

Measurements  of  growth  velocity  fall  into  two  general  classes:  those  in  which  free  dendritic- 
growth  is  measured  and  those  which  confine  the  growing  ice  in  a  tube  that  is  placed  in  a  bath  of 
some  other  liquid.  Both  methods  suffer  from  a  disadvantage  because  the  interface  temperature 
can’t  be  directly  measured. 

Measurements  of  the  growth  velocity  in  the  direction  normal  to  the  basal  plane  were  made  by 
Ifillig  (1958)  and  Sperry  (1965).  In  both  cases  the  growing  ice  was  confined  to  a  capillary  tube. 
Although  the  results  are  not  in  precise  agreement,  they  suggest  that  ice  grows  in  the  e-axis  direction 
by  surface  nucleation.  Damaged  crystals  grew  faster  than  undamaged  crystals;  they  thought  that 
damage  to  the  interface  provided  the  imperfections  necessary  for  dislocation-type  growth.  They 
found  that  a  minimum  supercooling  level  was  necessary  before  growth  was  observed -additional 
evidence  of  a  surface  nucleation  mechanism. 

Measurements  ot  the  growth  velocity  in  the  direction  parallel  to  the  basal  plane  have  been  made 
by  a  large  number  of  investigators  (Fletcher  1970.  Hobbs  1674.  Kallungal  and  Barduhn  1977). 
Fletcher  describes  the  results  of  many  investigators  as  showing  rough  qualitative  agreement.  All  the 
measurements  can  be  described  by  a  relation  of  the  form 

Gk  =  a(0j)6.  CD 

Unfortunately,  the  value  of  the  exponent  b  has  been  found  to  vary  front  1 .3  to  2.2,  and  measure¬ 
ments  made  at  the  same  value  of  0i  vary  widely. 

It  is  difficult  to  assess  how  accurately  the  data  represent  the  interface  kinetics  and  not  complica¬ 
tions  of  the  heat  transfer  effects.  The  results  were  generally  interpreted  as  implying  a  dislocation 
mechanism,  but  this  interpretation  has  been  criticized  by  Jackson  et  al.  (1967).  Kallungal  and 
Barduhn  (1977)  state:  “Interpreting  growth  rates  in  capillaries  is  a  difficult  and  unrewarding  task 
since  the  rates  are  as  much  dependent  on  the  properties  and  physical  dimensions  ol  the  capillary 
tube  as  the  characteristics  of  ice  and  water.”  To  resolve  the  question  ofa-axis  growth  rate,  they 
made  a  number  of  measurements  of  free  dendritic  growth  in  both  flowing  and  quiescent  water. 

They  found  that  at  large  impressed  water  velocities,  heat  transfer  of  the  boundary  layer  type  con¬ 
trolled  the  growth.  No  influence  of  interface  kinetics  was  found  to  exist  up  to  impressed  velocities 
of  68  ctn/s.  They  demonstrated  that  in  quiescent  water  the  growth  direction  with  respect  to  giavitv 
strongly  affects  the  results.  This  they  interpreted  as  evidence  of  heat  transfer  complications  due  to 
natural  convection.  The  importance  of  natural  convection  to  ice  growth  rates  was  also  demonstrat¬ 
ed  by  Gilpin  ( 1976).  The  rest  of  the  literature  contains  no  consideration  of  natural  convection 
heat  transfer. 
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In  summary,  the  different  growth  tales  parallel  to  the  basal  plane  and  normal  to  the  bus.ii  plane 
of  an  ice  crystal  are  the  result  of  differing  interlace  kinetics,  drouth  in  the  casts  piobahU  pioceed- 
by  surface  nucleation  for  perfect  crystals  and  by  a  dislocation  mechanism  tor  damaged  cixo.ilv  I  he 
interface  kinetics  of  (/-axis  growth  has  not  been  completely  defined:  the  mechanism  is  piobablv  that 
of  continuous  growth.  However,  it  appears  that  the  kinetics  are  very  fast  and  that  lot  practical  pm- 
poses  the  growth  rate  is  totally  controlled  by  the  rate  at  which  heat  is  transported  away  liom  the 
interface. 


Combined  effects  of  kinetics  and  heat  transfer  on  the  growth  rates 

The  rate  of  growth  of  an  ice  crystal  depends  on  two  piocesses  that  take  place  in  senes  the  at 
tachment  and  rearrangement  of  water  molecules  on  the  crystal  surface  in  accoidance  with  i  he  in¬ 
trinsic  kinetics,  and  the  transport  of  the  latent  heat  of  fusion  away  from  the  crystal  surface.  De¬ 
pending  on  the  relative  rates  of  these  two  processes,  either  could  control  the  glow  th  rate.  Prediction 
of  the  growth  rate  involves  determining  how  the  separate  processes  interact  and  how  the  two  can  be 
combined  to  give  an  overall  growth  rate. 

The  growth  rate  will  be  modeled  here  as  two  processes  m  series,  and  the  influence  o!  the  surface 
curvature  will  also  be  included.  Let  7.  be  the  bulk  temperature  of  the  supercooled  liquid,  and  7nl 
be  the  equilibrium  temperature  of  the  ice-water  mixture.  For  pure  water,  the  equilibrium  tempera¬ 
ture  will  be  0°C.  The  overall  supercooling  level  0  can  be  obtained  as  the  sum  of  the  three  tempera¬ 
ture  differences  representing  the  driving  forces  required  to  overcome  the  surface  curvature,  intrinsic 
crystallization  and  the  heat  transfer  resistances: 


°  =  Tm-'Tf=(Tm-T<:)nTe-Ti)UTrTf) 


where  (7'm-7'1)  is  the  temperature  difference  required  to  overcome  the  surface  curvature  resistance, 
(Te-Tj)  is  required  to  overcome  intrinsic  kinetics  resistance  and  (  T-  Tf)  is  required  to  overcome  the 
heat  transfer  resistances. 

The  surface  curvature  resistance  is  a  function  of  the  radius  of  curvature  of  the  ice/water  interface 
The  temperature  needed  to  overcome  this  resistance  can  be  estimated 


T  =  T 
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where  y  =  ice/water  interfacial  tension  (7  ^  22  erg  cm"2  [Fletcher  D>70] ) 

/.  =  latent  heat  of  fusion  per  gram  of  ice 
r  =  radius  of  curvature 
Pj  =  density  of  ice. 

Te  can  now  be  estimated  as  3  .Ox  I  O’6  cm  °('/ r.  Therefore,  if  r  >  ~  10" 1  cm  ( 1 0  pm)  then  the  curva¬ 
ture  effect  will  be  essentially  negligible.  We  can  assume  that  for  the  purposes  of  this  study  that  this 
will  always  be  the  case.  Therefore,  eq  22  can  be  rewritten  as 


0=7m-7"f=(7nl-7-i)  +  <7'j-7»  (24) 

where  (Tm-Tx)  is  required  to  overcome  the  intrinsic  kinetics  resistance.  Let  </  be  the  rate  of  heat 
transfer  per  unit  surface  area.  Assuming  a  steady  growth  rate,  we  see  that 

<7  =  <«Tm-r./  =  f,(Tr  7f)  =  /T(rnl- Tf)  (25) 


where  a  and  h  are  the  crystallization  kinetics  coefficients,  h  is  the  heat  transfer  coefficient  and  h 
represents  an  overall  transfer  rate.  We  can  make  the  crystallization  kinetics  linear  as 
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The  growth  rate  along  the  a-axis  determines  the  major  linear  dimension  of  the  fra/il  dish.  There¬ 
fore,  it  is  thea-a.xis  growth  rate  that  appears  in  the  number  continuity  equation  is  the  function  U. 
Kallungal  and  Barduhn  ( 1 977)  measured  zr-axis  growth  rates  and  found  no  evidence  of  a  rate  limit¬ 
ing  kinetic  step.  This  implies  that  a"  >  h.  Thus 
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Tire  growth  rate  in  the  u-axis  is  strictly  controlled  hy  heat  transfer.  An  expression  for  the  heat 
transfer  coefficient  h  will  be  determined  in  the  next  section. 

Growth  along  the  c-axis  is  much  slower  than  <r-axis  growth  for  all  si/es  of  crystals.  This  implies 
that  c-axis  growth  is  controlled  by  the  intrinsic  kinetics  and  thus  h  >a".  The  c-axis  growth  rate 
does  not  appear  in  the  number  continuity  equation.  The  latent  heat  released  by  growth  along  the 
c-axis  may  contribute  somewhat  to  the  overall  heat  balance:  however,  research  has  revealed  that  the 
latent  heat  released  by  c-axis  growth  is  effectively  negligible.  Therefore,  only  the  growth  along  the 
iz-axis  will  be  considered. 


Heat  transfer  from  ice  crystals  suspended  in  turbulent  water 

Introduction 

In  this  section  expressions  for  the  rate  of  heat  transfer  from  suspended  ice  crystals  will  be  formu¬ 
lated.  To  determine  the  transfer  ratio,  it  is  necessary  to  describe  the  ambient  velocity  distributions 
of  the  fluid  about  the  crystal.  Frazil  is  created  and  develops  only  in  water  that  is  turbulent.  Rivers 
and  channels  are  inherently  turbulent  because  of  the  instability  of  their  bulk  currents.  Wind  can 
make  large  water  bodies  become  turbulent.  Frazil  is  also  created  in  crystallizers  in  which  the  water 
is  made  turbulent  by  impellers,  turbines  or  other  means.  To  describe  the  velocity  distribution  of 
the  water  surrounding  the  crystals  requires  knowledge  of  the  properties  and  characteristics  of  turbu¬ 
lence.  The  first  part  of  this  section  will  be  a  very  brief  review  of  the  Kolmogorov  theory  o!  locally 
isotropic  turbulence.  For  further  details  the  reader  should  examine  the  texts  of  Batchelor  ( I '>53 ). 
Hinze  ( 1 959)  and  Tennekes  and  Lumley  ( 1972). 

Analytical  expressions  for  heat  transfer  have  been  developed  for  particles  immersed  in  a  station¬ 
ary  tluid.  in  a  fluid  moving  with  a  uniform  translational  motion,  and  in  a  tluid  whose  velocity  varies 
linearly  with  position  (shear  flow).  These  analytical  expressions  will  be  discussed  later  in  this  section 

If  the  ambient  velocity  about  a  suspended  particle  in  turbulent  water  cannot  he  described  in 
terms  of  these  velocity  distributions,  an  analytical  expression  may  not  be  possible .  In  this  case  a 
more  empirical  expression  is  necessary.  Frazil  crystals  are  subject  to  gravitational  and  meitial 
forces  which  give  them  an  additional  translational  motion  relative  to  tne  tluid  because  the  density 
of  ice  is  different  from  that  of  water.  The  magnitude  of  this  translational  motion  is  determined  in 
this  section  and  its  possible  influence  on  the  transfer  rates  is  assessed  Also,  as  the  ice  ciystals  me 
not  spherical,  the  influence  of  their  disk  shape  on  the  transtcr  rates  is  deteimmed.  Foi  the  growth 
of  frazil  in  fresh  water,  only  the  transfer  of  heat  must  be  determined.  The  corresponding  treatment 
and  results  for  mass  transfer  from  suspended  particles  are  identical  when  expiessed  appiopnatelv . 
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Fluid  turbulence 

The  physical  basis  ot  the  Kolmogorov  theory  of  turbulence  can  be  visualized  as  numerous  inter¬ 
acting  eddies  ot  all  possible  scales.  The  very  largest  eddies  originate  directlv  fiom  the  instabilities 
of  the  mean  bulk  How.  The  scale  and  orientation  of  these  largest  eddies  are  imposed  by  the  geome¬ 
try  ot  the  flow  situation.  In  a  stirred  crystal  I  i/ci  the  largest  eddies  are  created  bv  the  impeller  and 
are  in  scale  with  the  impeller  width  these  large  eddies  constitute  the  bulk  flow.  In  a  river  or  chan¬ 
nel  the  size  ot  the  largest  eddies  are  limited  by  the  depth  or  w  idth  of  the  channel.  In  a  large  body 
ot  still  water,  in  which  the  turbulence  is  generated  by  the  shear  stress  ot'  the  wind,  n  is  not  possible 
to  easily  predict  the  scale  ol  the  large  eddies,  unless  some  nonhomogeneous  feature,  such  as  a  iher- 
mocline.  exists. 

Knergy  is  extracted  from  the  large  eddies  through  the  inertial  interaction  of  these  eddies  w  ith 
smaller  eddies.  The  amount  ot  kinetic  energy  per  unit  mass  in  the  large-scale  eddies  is  proportional 
to  u~ ,  where  it  is  their  velocity.  This  energy  is  assumed  to  be  lost  in  a  time  pioportional  to  y  u. 
where  C  is  the  length  scale  ol  the  eddies.  Thus,  the  rate  at  which  energy  is  supplied  from  the  large- 
scale  eddies  to  the  smaller  eddies  is  u2u/i  =  t/’  ./t.  This  energy  cascades  through  the  spectrum  ol' 
eddy  sizes  to  the  smallest  eddies.  As  the  eddy  size  becomes  smaller,  the  geometric  orientation  of 
the  large  eddies  is  lost.  The  turbulence  is  isotropic  when  smaller  eddies  are  randomly  orientated. 
Tire  energy  cascade  is  not  affected  by  the  fluid  viscosity  until  the  smallest  scales  are  reached,  where 
this  energy  is  dissipated  by  the  viscosity.  Tire  dissipation  rate  must  equal  the  rate  at  which  energy 
is  supplied  to  the  small-scale  eddies.  Therefore,  the  dissipation  rate  e  can  be  defined  as 


e  =  ir'/V.  (29) 

The  scale  at  which  viscous  dissipation  becomes  important  can  be  estimated  if  the  fluid  viscosity  v 
and  e  arc  known.  From  these  parameters  it  is  possible  to  form  a  length  scale  r}  such  that 

r]  ~~(v3/e)'li  .  (30) 

r?  is  the  dissipation  length  scale  or  the  Kolmogorov  scale. 

For  Hows  with  a  sufficiently  high  Reynolds  number,  if  the  spectrum  of  eddy  sizes  is  normalized 
by  the  dissipation  length  scale,  there  is  a  universal  character  for  the  range  of  eddy  sizes  smaller  than 
the  large  energy -containing  eddies.  This  is  the  universal  equilibrium  range.  The  universal  character 
of  many  different  turbulent  flows  is  shown  in  Figure  4.  which  is  plotted  to  show  the  energy  of 
eddies  against  the  wave  number  of  the  eddies  (reciprocal  of  wave  length  ). 
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3.8x104  grid  turbulence 
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Figure  4.  Universal  turbulence  spectra  (after  Wutiia  I  Vs). 
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The  range  of  scales  larger  than  the  dissipation  scale  hut  miuIIci  than  the  wale  ■  >!  the  huge  cnctex 
containing  eddies  is  the  inertial  subrange.  It  is  impossible  to  picdict  a  pnoii  ibe  maximum  se.de  al 
which  the  inertial  subrange  will  begin.  Within  the  menial  suhiange  the  blind  uscnsiu  lias  no  elicn 
and  all  energy  dissipation  results  from  the  inertial  interactions  between  eddies  ol  dilleient  si/es 
Therefore,  the  only  scaling  parameter  available  is  t .  which  can  be  interpreted  as  the  tale  (pei  unit 
mass  ol  fluid)  at  which  energy  cascades  through  the  spectrum  of  eddies  si/es  within  the  inertial 
subrange. 

The  range  of  scale  that  is  smaller  than  the  dissipation  scale  is  the  dissipation  suhiange.  The  fluid 
viscosity  play  s  an  important  role  in  the  subrange  and  acts  quickly  to  dampen  and  dissipate  the 
fluid  motion.  This  small-scale  motion  automatically  adjusts  itself  to  the  value  of  the  viscosity  and 
the  rate  of  energy  transfer. 

Heat  transfer  from  suspended  particles 

This  discussion  of  heat  transfer  from  force-free  particles  will  concentrate  on  two  kinds  ol’  flows: 
first,  particles  suspended  in  fluid  with  steady  velocity  distribution:  second,  particles  suspended  in 
turbulent  fluid. 

Linear  velocity  distributions.  Probably  the  most  basic  situation  that  can  be  anaiv/ed  is  that  of 
an  isothermal  particle  suspended  in  a  motionless  fluid,  with  no  relative  motion  between  the  particle 
and  the  fluid.  In  this  case  the  heat  transfer  rate  is  controlled  purely  by  the  rate  at  which  beat  can 
diffuse  from  the  surface  of  the  particle  to  the  bulk  of  the  fluid  (Carslaw  and  Jaegar  I ‘>59).  What¬ 
ever  the  particle  shape,  the  steady  distribution  of  temperature  due  to  diffusion  becomes  spherically 
symmetric  at  large  distances  from  the  particle.  At  these  distances  the  temperature  distribution  is 
the  same  as  that  which  would  be  caused  by  continuous  point  sources  of  heat,  emitting  beat  at  the 
same  rate  as  the  actual  particle. 

It  is  convenient  to  define  the  nondimensional  measure  of  the  heat  transfer  rate  or  Nusselt  num¬ 
ber  Nu  as 

Nu  =  (hr  Ik)  (51) 

where  r  is  the  major  radius  of  the  particle,  k  is  the  thermal  conductivity  of  the  fluid,  and  h  is  the 
heat  transfer  coefficient.  For  a  spherical  particle  of  radius  /  in  a  stationary  fluid,  the  heat  transfer 
is 

4nr  k(Ts-  T()  =  4nr2  h(Ts-Tf)  (32) 

where  Ts  is  the  temperature  at  the  particle  surface  arid  Tf  is  the  bulk  temperature  of  the  fluid. 
Therefore 

Nu0  =  I  •  (33) 

The  method  used  to  determine  the  transport  rate  from  the  particle  suspended  in  a  fluid  moving 
with  a  steady  uniform  velocity  will  depend  on  the  value  of  the  Peclet  number,  defined  as 

Pe=(Vfr/a)  (54) 

where  Vf  is  the  steady  translational  velocity  and  a  is  the  thermal  diffusivity  of  the  fluid. 

If  the  Peclet  number  is  small  (Pe  <  I ),  the  transport  from  the  particle  surface  is  dominated  by 
diffusion.  The  effect  of  the  fluid  motion  is  to  modify  the  spherically  symmetrical,  steady-state 
distribution  of  temperature  due  to  diffusion  at  large  distances  from  (lie  particles  and  to  increase 
the  transport  rate.  Because  tl*e  influence  of  the  fluid  motion  is  only  evident  at  large  distances  liom 
the  particle,  the  transport  rate  is  relatively  insensitive  to  the  Reynolds  number  and  the  particular 
form  of  the  ambient  flow  field. 


The  solution  of  the  steady-state  heat  transport  equation  I'rom  a  suspended  paitiele  at  a  small 
Peelet  number  belongs  to  the  class  of  problems  known  as  singular  per  turbnlion  pmblemv  I  lie 
steady-state  problem  was  solved  by  Aerivos  and  Taylor  (  |9<>2 )  and  Brcnnei  ( I  l)o  ’>  |  |n  using  "innei  " 
and  “outer”  expansions  of  the  temperature  fields  and  matching  these  solutions  m  then  common 
domain  of  validity.  The  inner  solution  held  for  the  region  near  the  particle  in  which  ditlusion  was 
assumed  to  dominate.  The  outer  solution  held  tor  the  region  tar  Imni  the  panicle  where  comes 
tion  was  assumed  to  dominate.  The  solution  for  the  transfer  rate  is 


NU-NUo  I  n. 

-NiT-2  N“° 


<3'l 


If  the  Peelet  number  is  large  (Pe  >  1 ),  the  transport  from  the  particle  surface  is  dominated  In 
convection.  The  gradients  of  temperature  exist  only  in  a  small  boundaiv  layer  near  the  parti Je. 
outside  of  which  the  temperature  can  be  assumed  to  be  uniform.  The  transfer  late  at  a  large  Peelet 
number  is  dependent  on  the  form  of  the  fluid  motion  near  the  particle  and  therefore  the  shape  ot 
the  particle.  A  uniform  translational  motion  that  is  linear  with  position  corresponds  to  the  condi¬ 
tions  of  Stokes  How.  The  expression  for  transfer  from  spherical  particles  in  Stokes  llow  at  large 
Peelet  numbers  has  been  determined  by  Levich  ( 1962).  Brian  and  Hales  ( 1909)  and  Balcheloi  I  1979). 
The  solution  of  Batchelor  is 


Nu  =  Nu0  +0.624  Pe1 


The  third  and  most  general  situation  is  that  of  a  particle  immersed  in  a  fluid  w  ith  a  steads  velo¬ 
city  distribution  that  varies  linearly  with  position.  This  type  of  llow  shear  (low  has  not  had  the 
extensive  analysis  of  the  previous  two  flow's.  In  general,  a  linear  velocity  distribution  can  be  repre¬ 
sented  as 


1  f  5ij  xi 


where  is  position  and  the  tensor  Sy  can  be  written 


Sjj  =£  +  S2 


(37) 


(38) 


where  E  represents  the  pure  straining  motion  ot  the  fluid  and  52  represents  a  solid  body  rotation  ot 
the  fluid.  The  angular  velocity  of  the  solid  body  rotation  is  determined  by  the  vorticity  ol  the  am¬ 
bient  flow.  We  can  see  that  a  5  of  any  magnitude  can  be  produced  by  an  infinite  number  of  com¬ 
binations  of  E  and  52. 

Which  combinations  of  E  and  52  are  relevant  tor  determining  the  transfer  rates  trom  a  suspended 
particle?  Batchelor  (1971)  provided  the  answer.  An  axis  of  symmetry  along  the  particle  can  be  de¬ 
termined  by  resolving  the  components  of  vorticity  along  the  principal  axes  ol  the  rate  ot  strain  ten¬ 
sor  E .  It  is  the  strain  rate  Ew  along  this  axis  of  symmetry  that  is  largely  responsible  lor  determin¬ 
ing  the  transfer  rates  from  a  particle.  Transfer  in  all  other  directions  will  be  suppressed  by  the  rota¬ 
tion  of  the  fluid. 

Tor  shear  flows  the  Peelet  number  can  be  defined  as 


Pe 


r2S 


( 3°  I 


a  a 

For  small  Peelet  n umbers,  Batchelor  ( 1 979)  determined  the  transfer  rate,  to  the  first  oidei .  as 


Nu-Nuo 


Nu0 


=  0.40  Nu0  Pe  . 


For  large  Peelet  numbers.  Batchelor  ( 1979)  determined  the  transfer  rate,  to  first  order,  as 


Nu  =  Nuo+0.97  Pe1 


(40) 


(41  ) 
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Vi '  nhnear  vch  ><.  // »•  Jistrtbuti .  >n  Pat  i  teles  suspended  hi  .i  ilniii  moving  w  tilt  .1  in  mi  In  k-.i  i  \cl.  n  in 
distribution  (relative  to  the  particle)  are  suh|cct  in  .ill  situations  in  wlikh  the  Revnolds  numbei  is 
large  or  the  fluid  is  turbulent  At  a  large  Rev  Molds  minihei  the  link!  houtulaiv  l.nei  Imins.  wjki 
separation  occuis.  the  houmi.nv  la\er  becomes  l  nr  t'nlen  i .  etc  1  Ins  n  impoi  l.iiit  tomato  iinlusin.il 
processes  because  tliev  depend  on  the  transfer  ot  heat  to  ot  tioin  an  object  aionnd  which  .1  lluid  n 
flowing.  It  the  shape  ot  the  obiect  is  complex,  the  tianslei  tales  imisi  be  detemniied  enipii  tcallv  . 

Il  the  shape  ot  the  object  is  simple,  such  as  a  single  c\  liiulei  01  spheie.  the  heat  tianslei  late  can 
be  determined  front  known  conelations  ot  heat  tianslei  tales  and  the  piopemo'  ot  the  flow  I  he 
Frossling  equation  is  a  well  known  correlation  that  relates  the  heat  tianslei  tioin  a  spheie  to  the 
Reynolds  number  and  Prandfl  ntimbei  of  a  steads  flow  Flits  telaiion.  which  has  been  loimd  to  be 
accurate  for  Reynolds  numbers  betw  een  I  and  10'.  n 

Nu  =  Nu0+0.4:  Re1  I’r1  (4_'| 

Re  is  the  Reynolds  number  defined  as  Re  -  r  I ,  t  Hie  I’landtl  numbei  is  dclmcd  as  t  o  (kinematic 
viscosity/thernial  diffusivity ). 

At  this  point,  expressions  tor  the  transtei  rate  Itom  panicles  suspended  111  a  tuibttlent  fluid  will 
be  determined.  In  a  previous  section  the  dissipative  and  inertial  subiatiges  ot  the  tiubulent  spec- 
truni  were  described.  We  will  see  that  ’he  ambient  veloctts  dtstiihution  ot  the  dissipative  subtauge 
can  be  described  in  terms  of  linear  motion  and  the  inertial  subrange  cannot ;  sery  d  1  tie  lent  means 
of  determining  the  transfer  rates  are  required  in  each  subrange  Thetefote.  these  two  subranges 
can  be  said  to  compose  two  regimes  of  heat  transfer  l  nlortunatelv  .  we  cannot  know  beforehand 
which  regime  the  t'ra/.il  will  be  in.  as  the  st/e  ot  the  crystals  and  level  ot  turbulence  max  vary  widely 
in  natural  water  bodies. 

Particles  in  turbulence  dissipative  regime.  It  the  avstal  st/e  is  small  relative  to  the  Kolmogorov 
length  scale,  it  is  in  the  dissipative  regime.  In  the  dissipative  regime  the  fluid  eddies  ate  strongly 
dampened  and  dissipated  by  the  fluid  viscosity .  In  effect .  the  crystal  is  smaller  than  the  smallest 
scales  ol  the  turbulent  eddies.  It  does  not  experience  the  turbulence  as  interacting  eddies  but 
rather  as  a  fluid  motion  that  varies  linearly  with  position.  The  magnitude  of  .S'  in  the  dissipative 
regime  can  be  estimated  as  (Tennekes  and  Luntlev  llP2) 

.V~(e/u)i:.  (4  s) 

Tire  Reynolds  number  ot  the  particle  motion  will  be  (assuming  the  particles  experience  ottlv  this 
shear) 

Re  =  (r2  c 1  2  )/i>'  2  (44) 

and  the  Peclet  number 

Pe  =  ( r 2  e1  2  )/(otv'  2 ).  145  ) 

As  r  <  r?  in  the  dissipative  region . 

r  <  (v'/c)'  4 
and  therefore 

Re  <  I. 

For  the  Peclet  number  to  equal  unite,  that  is.  Pe  =  I. 


r  =  (1/Pr1  :)r? 


(4(i ) 


where  Pr  is  the  Prandtl  number.  Therefore,  when 
r  <  ( 1  / Pr 1  2  )rj:  Pe  <  1 

r  >  (l  /Pr1  2  )rj;  Pe  >  1 . 


(47a) 

(47b) 


In  the  analysis  of  mass  transfer,  usually  the  Praiultl  number  (or  analogously  for  mass  transfer, 
the  Schmidt  number)  is  much  larger  than  1 .  Therefore,  the  low  Pc'clet  case  corresponds  only  to 
very  small  particles,  usually  much  smaller  than  the  range  of  interest.  However,  for  heat  transfer 
from  frazil  crystals,  the  Prandtl  number  is  approximately  I  3.  and  the  low  Peclet  case  is  relevant. 

It  is  interesting  to  note  that  there  is  only  a  narrow  range  in  which  Pe  >  I  and  Re  <  I .  As  the  parti¬ 
cle  size  becomes  large,  r  will  approach  the  Kolmogorov  scale,  rj.  and  the  conditions  of  the  dissipa¬ 
tive  regime  will  no  longer  be  valid.  However,  data  suggest  (Batchelor  1980)  that  the  tlow  distribu¬ 
tion  can  be  considered  linear  until  r  ==  I  Op.  which  expands  the  large  Peclet  number  range.  In  the 
interest  of  generality,  then,  the  result  for  both  the  small  and  large  Peclet  cases  will  be  given. 

As  noted,  Batchelor  identified  /.w  as  the  important  component  of  the  shear  for  determining  the 
transfer  rates.  The  average  over  time  <|/fw|>  is  a  parameter  of  the  turbulent  fluid  in  which  the 
particle  is  immersed  and  is  independent  of  the  properties  of  the  particle.  The  small-scale  properties 
of  the  turbulence,  in  particular,  determine  its  value. 

Batchelor  (1080)  showed  that  in  locally  homogeneous  and  isotropic  turbulence 


<|£wl>  =  0.18(e/u)'  2. 


(48) 


Therefore,  Batchelor’s  result  for  the  heat  transfer  at  a  small  Peclet  number  (eq  40)  can  be  written 


Nu-Nup 

Nu0 


=  0.1 7  Nu, 


lrl 

'U  u12/  • 


(49) 


The  large  Pe'clet  number  result  of  Batchelor  (eq  41 )  can  be  written 
/  r 2  6 1  /2  \ 1  2 

Nu  =  Nuo+0.55(— Tj-j  .  (SO) 


These  results  are  for  a  particle  immersed  in  a  shear  flow  produced  by  the  turbulence  of  the  fluid. 

It  is  assumed  that  there  are  no  inertial  forces  or  buoyancy  forces  acting  on  the  particle  that  would 
cause  additional  movement  of  it  relative  to  the  fluid. 

Particles  in  turbulence- inertial  regime.  If  the  crystal  size  is  large  relative  to  the  Kolmogorov 
length  scale  it  is  in  the  inertial  regime.  A  number  of  theories  of  mass  transfer  in  the  inertial  regime 
have  been  developed  and  abandoned  by  the  chemical  engineers.  These  are  discussed  by  Wadia 
(1975).  The  cause  of  the  difficulty  is  that  the  flow  field  in  the  vicinity  of  the  crystal  is  complex 
and  the  transfer  may  proceed  at  many  scales.  The  particle  may  interact  with  fluid  eddies  that  are 
both  larger  and  smaller  than  it  is.  The  problem  then  becomes  characterizing  the  ambient  velocity 
distribution  so  as  to  accurately  determine  the  heat  transfer  rate. 

The  ambient  velocity  can  be  characterized  in  many  different  wavs,  each  corresponding  to  a 
different  eddy  size.  It  seems  reasonable  to  assume,  following  Wadia  ( 1975).  that  the  predominant 
shear  that  the  particle  will  experience  will  be  produced  by  eddies  closest  to  the  particle  .that  are  of 
the  same  size  as  the  particle.  Hddies  that  are  significantly  larger  than  the  particle  will  entrain  both 
the  particle  and  the  fluid  around  it.  Very  small  eddies  relative  to  the  particle  size  may  enhance  the 
overall  transport  by  some  mechanism  of  renewal  of  the  boundary  layer  surrounding  the  particles. 

IT 


but  it  is  eddies  of  a  size  comparable  to  the  si/e  of  the  particle  that  will  cause  the  most  significant 
gradients  near  the  crystal  surface.  Therefore,  the  shear  can  be  estimated  as  (Levich  1962) 


S~elir2i.  (51 

This  shear  will  not  be  linear  with  position.  The  Reynolds  number  of  the  particle  motion  will  be 


and  the  Peclet  number 


As  r  >  17  in  the  inertial  regime 
r  >  (ir'/e)'  4 
and  therefore 
Re  >  1 

and  the  Peclet  number  is 
Pe  >  Pr. 

In  the  inertial  subrange,  a  linear  ambient  velocity  distribution  does  not  exist.  The  Peclet  number 
is  large,  therefore  gradients  of  temperature  will  exist  only  in  small  boundary  layers  near  the  particles. 
Fluid  boundary  layers  will  also  exist  and,  if  the  Reynolds  number  becomes  large  enough,  they  may 
become  turbulent.  This  situation  cannot  be  analyzed  in  terms  of  the  linear  velocity  distributions 
described  earlier.  Therefore,  to  determine  the  heat  transfer  from  the  particles,  the  Frossling  equa¬ 
tion  will  be  used,  following  the  example  of  Wadia  (1975).  The  mean  square  shear  between  two 
points  separated  by  distance  r  in  the  inertial  subrange  can  be  estimated  as  (Batchelor  1953) 

\/¥  =  2J  (e'3)r2  3.  (54) 

Substituting  eq  54  and  52  into  eq  42,  we  see  that  the  Frossling  equation  becomes 


.  =  Nu0+0.7o(-  iy£---)  Pr'  3- 


The  applicability  of  this  equation  to  the  situation  under  discussion  here  remains  open  to  question. 
This  equation  was  derived  for  steady,  nonturbulent  flow,  not  shear  flow.  However,  it  is  interesting 
to  note  that  if  the  size  of  particle  r  is  made  into  a  number  without  dimensions  by  the  Kolmogorov 
scale  such  that  r*  =  r/77,  then  the  Frossling  equation  and  the  expression  for  the  transfer  rate  with  a 
large  Peclet  number  and  a  general  linear  ambient  velocity  distribution  reduce  to 

Nu  «  (r*)2  3  Pr13  (56) 


and  therefore  have  the  same  dependence  on  the  independent  parameters  of  size  and  the  Prandtl  number . 
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We  write  the  Frbssling  equation  using  the  eltaraetettstie  slieai  of  an  eddy  si/e  equal  to  the  pani¬ 
cle  si/e,  so  wltat  can  he  determined  about  the  influence  ol  eddies  larger  and  smaller  than  the  pani¬ 
cle  size'.’  Wadia  ( 1  L> 7 5 )  calculated  that  eddies  smaller  than  hall  the  panicle  diametei  account  lot 
less  than  20'-'  of  its  relative  velocity  .  lie  determined  that  eddies  largei  than  the  particle  diameivi 
accounted  for  657  of  its  relative  motion,  and  eddies  greater  than  live  particle  diameters  were  re¬ 
sponsible  for  less  than  1  2''i .  kuboi  et  al.  ( I  ‘■>74 )  used  high  speed  photographs  to  study  the  motion 
of  neutrally  buoyant  trace  particles  suspended  in  a  turbulent  stirred  tank  and  in  a  flowing  pipe. 
They  found  that  the  mis  velocities  of  the  particles  were  2.0  ( t  r ) 1  ' .  Tins  supports  the  theory  that 
it  is  the  eddies  of  si/e  comparable  to  that  of  the  particle  that  cause  the  Hind  motion  around  the 
particle. 

The  small-scale  motion,  smaller  than  the  si/e  of  the  crystal,  may  enhance  the  heat  transput) 
from  the  crystal  by  penetrating  the  boundary  layer  around  the  crystal.  It  is  difficult  to  quantify 
this  process  but  this  enhancement  has  been  successfully  accounted  for  (although  empirically)  by 
correlation  of  the  turbulent  intensity,  a,  .  of  the  fluid,  a,  is  defined  as 

a  T  =  \ JiC  I  lr  (57) 

where  \Jir  is  the  mis  value  of  the  velocity  deviation  from  the  mean  velocity  l'{.  The  experiments 
of  Lavender  and  Pei  ( 1%7)  demonstrated  that  the  I  rdssling  equation  could  be  written  as 

Nu  =  Nuo+0.44aal  Re<,/2  >a)  Pr1  3  (58) 

where  a  is  an  experimentally  derived  coefficient.  They  found  that  for 

aT  Re  <  1000.  a  =  0.045 
and 

aT  Re  >  1000.  a  *=  0.25 

with  the  break  occurring  at  the  point  at  which  the  boundary  layer  becomes  turbulent.  Therefore, 
having  aT  Re  <  1 000  gives 

Nu  =  Nu0  +070a°  035  Re0535  Pr1/3  (60a) 

and  having  aT  Re  >  1 000  gives 

Nu  =  Nu0  +0.70a°-25  Re°-7S  Pr1/3.  (60b) 

For  a.  Re  >  1000.  Nu  is  essentially  independent  of  the  crystal  size.  This  range  is  where  industrial 
crystallizers  generally  operate  and  the  independence  of  the  transfer  processes  and  particle  size  is 
often  seen.  This  observation  has  been  generalized  and  is  called  McCabe's  Law  (see  the  Steady -Stau 
Crystal  Number  Continuity  Equation  and  Heat  Balance  in  a  MSMPR  Crystallizer  section ).  This 
correspondence  between  theory  and  observation  is  encouraging. 

Heat  transfer  from  suspended  ice  crystals 

At  this  point  the  heat  transfer  equations  developed  in  the  preceding  section  will  be  applied  to 
ice  crystals  suspended  in  turbulent  water.  As  the  density  of  ice  crystals  is  different  from  that  of 
water  [(Pj-pf/pf)  =  87  | ,  the  crystals  are  subject  to  gravitational  and  inertial  forces  that  give  them 
a  translational  motion  relative  to  the  fluid.  The  magnitude  of  this  translational  motion  must  be 


determined  and  it;,  possible  Influence  on  the  transfer  tales  assessed.  Also,  as  the  ice  ci >  stals  arc 
not  spherical,  the  influence  of  then  disk  shape  on  the  transfer  rates  must  be  detei mined. 

Relative  translational  motion.  Let  the  translational  motion  of  a  ctystal  relative  to  the  fluid  be 
(  K  .  It  rS  l  -  I .  where  .S'  is  the  shear  of  the  ambient  tluid.  then  the  transfer  rates  are  detemimed 
pnmartK  on  the  basts  of  the  relative  translational  motion.  The  influence  of  the  shear  would  be  a 
second  order  or  higher  effect.  If  rS'l  (<  P  I .  then  the  transfer  rates  are  determined  primarily  on 
the  basis  of  the  shear  of  the  ambient  fluid.  When  rS  t  R  *  I .  Wadia  ( l‘>75 )  assumed  that  an  ef  fec¬ 
tive  crystal  velocity.  I  K  eff,such  that 

rKcff  -  K r +  r2R  1 1  : 


could  be  used  to  calculate  the  transfer  rate.  The  work  of  Batchelor  ( 1 1>7S).  Ih80)  supports  this 
means  of  combining  the  velocity  due  to  shear  and  that  due  to  translational  motion  in  estimating 
transfer  rates. 

Inertial  forces.  The  relative  motion  due  to  the  inertia  of  ice  crystals  will  now  be  determined. 
A  force  balance  on  a  crystal  can  be  written  as 


Jle  dVf  nt  Vr  dVo 

Wi  -jr  -  'n+V’Pr  ~jr 


dt 


(62) 


where 


A^  =  volumetric  shape  factor  (A'vr'*  =  volume  of  crystal) 
I  f  =  velocity  of  the  tluid 
V  =  velocity  of  the  crystal 
Pi  =  density  of  the  crystal 
p(  =  density  of  the  fluid 
X  =  added  mass. 


We  will  assume  that  1 )  gravity  effects  are  negligible  (they  will  be  analyzed  separately),  2)  the  density 
of  crystals  is  small  enough  so  that  the  tluid  properties  are  essentially  uninfluenced  by  the  presence 
of  the  crystals,  3)  all  the  elements  of  an  eddy  are  accelerated  identically,  4)  the  Basset  term  can  be 
ignored,  and  5)  the  tluid  turbulence  is  statistically  steady.  Now  let 


A{  =  (dVf/dr)  (63a) 

AH=[dVf/dt)-(dVJdt)  (63b) 


Vr 


(63  c) 


then 


Pi  AV r '  A  k  =(Pj'Pf)A'vr1z(f-/lK  Kvr3pfx-Fu  (b4) 

as  the  vectors  Af  and  AH  are  parallel. 

To  proceed  from  this  point,  the  size  of  the  crystal  relative  to  the  scale  of  the  turbulence  must 
be  known. 

If  the  crystal  is  in  the  dissipative  range,  then  r  <  r?  and  the  drag  force  can  be  estimated  using 
Stokes  Law 


C’„  =  12/Re. 


(65) 


where  Cr)  is  the  drag  coefficient.  Thus 


The  relative  acceleration  and  the  fluid  acceleration  can  be  estimated  in  the  following  manner 


-If 


u*,  r* 


((>7) 


where  T*  is  the  characteristic  time  of  the  eddy  of  interest  and  u*  the  characteristic  velocity.  As 
r  <r\  these  quantities  are  estimated  at  the  scale  of  17.  Thus 


(e>)' 


(68) 


and 


AR=(UR/T*)  =  UR(e/v)'  2. 

Substituting  into  eq  66  and  solving  for  UR ,  we  find 

( Kvr2/bnv )  |(prpf)/pf|  (e3/v)'  4 
Ur  (Kvr2/ 6irv)  [(Pj+pf X)/pf\  (e/u)12+l  ' 

Now  the  magnitude  of  the  relevant  shear  rate  can  be  estimated  as 

S  =  0.18  (e/u)1 2 

and  the  relative  magnitudes  of  UR  and  S  can  be  estimated  as 

rS  _  (0-18A~v/6ff)  [(p(+pf X)lpf |  (r/vf  + 1 
(A'v/6rr)  [(prp{)/pf)(r/v) 

and  it  can  be  seen  that  if  r  <  77 


(69) 


(70) 


(71) 


(72) 


£>■ 


(73) 


where  UR  is  the  relative  translational  motion  produced  by  the  inertial  forces  on  the  crystal. 

The  inertial  force  on  crystals  in  the  inertial  regime  will  now  be  determined  following  the  exam¬ 
ple  of  Levich  (1962)  and  Wadia  (1975).  Starting  with  eq  64,  we  can  estiniate  the  relative  accelera¬ 
tion  between  the  crystal  and  the  fluid  as 


ak~ur/t* 


(74) 


where 


t*  =  e/t/R 


(75) 


and  ?  is  the  characteristic  eddy  length  scale.  Thus 


To  characterize  the  velocity  of  the  fluid,  only  the  size  of  the  eddy  2  and  the  dissipation  rate  c  are 
available.  Therefore 


u*  =  (e2) 


(77) 


and 


A{*u*/T* 


(78) 


where 


T**i/u*. 


(79) 


Titus 


/lf  =  (w*2  T*)  =  elaV'a . 


(80) 


F0  can  be  written 


Fn=(C^I2)r2  pfU2R 


(81) 


then 


(prpf)*vr3(e2)2'3 

Ur  p,Kvr3  +  pfKvr3XHCu!2)nr2  pf2 


(82) 


We  can  see  that  (JR  is  a  function  of  C.  At  a  certain  value  of  2  =  2*.  the  relative  velocity  has  a 
maximum  value  Uu  .  2*  can  be  determined  such  that 

Kmax 

dUHl 32  =  0.  (83) 


2*  is  then 


*_  ^PjKy+PfK^X) 

?  (Cu/2)rrr2pf 

Substituting  into  eq  82,  we  find 

u  =j_/ _ _ y 

Kmax  dS  A\(pj+pfA')1  3(Cupf)2  V 


(er)' 


(84) 


(85) 


As  Levich  (1962)  found,  the  numerical  coefficient  of  eq  85  probably  has  no  great  significance. 
It  is  provided  to  indicate  the  absence  of  large  numerical  coefficients. 

We  recall  that  in  the  inertial  regime  rS  =  2.7(er)'  3;  thus 


rS  t(Pi+PfX)r3(Cnp1)2  Y 


■F 


(Pj-Pf) 

1 '  m  a  v 


and.as  1 2  <  Cn  <  1 .1  and  (pj-pf )  =  0.08,  we  can  see  that  when  r  >  rj 
(rS/UR  )  >  1 


(86) 


(87) 
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where  Lf  is  the  maximum  relative  translational  motion  pioduced  In  the  inertial  loreeson  the 

anj\ 

crystal. 

We  can  conclude,  therefore,  that  to  a  pood  approximation  the  inertial  force  can  be  ignored  in 
determining  the  transfer  rates  from  suspended  ice  crystals. 

Gravity  J forces.  The  relative  translational  motion  due  to  the  force  of  gravity  can  be  estimated  by 
equating  the  drag  force  and  the  gravity  force: 

P~K^r'^  =  \  rrr'Cnpfl- 
P  f 

where  A'v  is  the  volumetric  shape  factor,  l is  the  terminal  rise  velocity  of  a  crystal,  and  Cn  is  the 
drag  coefficient.  The  drag  coefficient  of  a  disk  whose  major  radius  is  perpendicular  to  the  (low  is 
a  well  known  function  of  the  Reynolds  number  and  is  available  in  many  texts.  Txperimental  ob¬ 
servation  has  shown  that  disks  under  the  influence  of  gravity  do  not  always  move  steadily.  At  the 
higher  Reynolds  numbers,  they  also  oscillate,  glide  and  tumble.  1’ach  of  these  motions  influences 
the  drag  coefficient.  Turbulence  also  influences  the  drag  coefficient  in  a  complex  manner.  In  gen¬ 
eral,  these  influences  will  all  tend  to  increase  the  drag  on  the  disks  slightly.  In  Figure  5  the  terminal 
rise  velocity  is  plotted  as  a  function  of  the  major  radius.  The  velocity  was  estimated  by  assuming 
that  the  disk  rises  steadily  with  its  axis  perpendicular  to  the  vertical.  This  may  be  an  overestima¬ 
tin'  of  the  rise  velocity. 

In  the  Stokes  range  of  r  <  0.03  cm 

Ut  =0M(g'v'  r2).  (89) 

In  the  intermediate  range  of  0.3  cm  OCO.14  cm 

L\  =  0.l6(g'°'715  u-0-428r'.i4)  (90) 

and  in  the  fully  turbulent  range  r  >  0.14  cm 

Ut  =  2'1  2  (g'r)'  2  (91) 

where  g’  is  the  reduced  gravity  (g-  2[(pf-Pj)/pf  |  [gK^/n] ).  To  compare  the  shear  motion  and  the 
relative  translational  motion  due  to  gravity,  several  assumptions  will  be  made  that  are  based  upon 
physical  knowledge  of  frazil  and  the  likely  sizes  of  the  frazil  crystals. 


Fully 
T  urbulent 
Range 


Intermediate 
I  Range 


Velocity  (cm/s) 

Figure  5.  Terminal  rise  velocity  of  frazil  disks. 


We  will  assume  that,  it  the  ice  crystal  is  in  the  Stokes  range,  that  is  r  <  0.0.'.  then  11  will  j|s<»  he 
in  the  dissipative  regime.  Then 

O'  -  Ut>)' : 

Ut  0.08^V 

and  using  the  best  estimate  for  v  and  we  see  that 


O  =  0.02  e1  3 
Ut  “  r 


r  <  0.02  cm. 


For  the  intermediate  range,  die  crystal  may  be  in  either  the  dissipative  or  inertial  regime.  In  the 
dissipative  regime 


rS  =  0.02  e1  2 

<7,  '  ,0.14  • 


0.02  cm  <  r  <  0. 1 4  cm 


and  in  the  in tertial  regime 


—  =  0.03  cm  <r< 0.14  cm.  (94) 

Ut  ro.8t 

The  fully  turbulent  range  corresponds  to  particles  with  radius  larger  than  0.14  cm.  This  is  very 
large  for  a  frazil  disk;  very  few  frazil  crystals  have  been  recorded  larger  than  this  size.  Therefore, 
this  range  will  not  be  considered  further. 

Given  the  maximum  sizes  of  disks  in  the  Stokes  range  and  the  intermediate  range,  e  need  only 
be  slightly  larger  than  1  for  the  shear  velocity  to  dominate.  This  value  of  1  corresponds  to  a  rela¬ 
tively  small  level  of  turbulence.  In  most  situations  it  will  be  possible  to  effectively  ignore  the  traits 
lational  motion  induced  by  gravity  in  determining  the  heat  transfer. 

Heat  transfer  from  disks.  Up  to  this  point,  the  heat  transfer  relationships  have  been  developed 
for  spheres.  The  influence  of  the  nonspherical  shape  of  ice  crystals  must  be  assessed. 

The  heat  transfer  from  a  disk  in  a  stagnant  fluid  has  been  determined  by  Wadia  ( 1975)  and 
many  others.  Wadia  demonstrated  that  if  the  major  dimension  of  the  disk  was  defined  as 
(^4 s/4 7r)'  7  =  r,  where  As  is  the  surface  area  of  the  disk,  then  for  a  disk  with  aspect  ratio  0.24, 

NuDIsk  =  <*»'/*)=  1.  (95) 

The  heat  transfer  from  the  edge  and  face  of  the  disk  can  be  represented  as 
Nuedck  =(hrjk)=  I 
Nu  hack  =(hrf/k)=  1 

where  re  is  one  half  the  thickness  of  the  disk  and  rf  is  the  radius  of  the  face  of  the  disk.  Noting 
that  r  =  0.92 r.  re  =  0.24r,  and  rf  -  r ,  we  find 

Nut )  ts  k  ~  (hr/k )  =  l.l 

Nufih;k  =  1  hr/k)  =  2.94 


Nui-  ack  =  (hr/k)  =  1 .0. 
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Figure  6.  Fdge  vs  face  transfer  (after  Wadia  19  751. 


Tiie  heat  transfer  rate  is  much  larger  for  the  edge  than  the  face  of  a  disk.  A  general  wa\  of  ex¬ 
pressing  this  is 


Nu  =  ( hr/k )  =  1 


(%l 


where  r  =  re,  r(,  r. 

If  the  crystals  are  in  a  turbulent  fluid  and  are  small  enough  such  that  r  <  pU’r)1  2 .  they  are  in 
the  small  Pe’clet  number  range  where  the  heat  transfer  is  dominated  by  diffusion.  At  small  Pe'clet 
numbers  the  heat  transfer  rate  is  insensitive  to  the  shape  of  the  particle,  as  noted  earlier.  There¬ 
fore  the  disk  shape  will  not  influence  the  heat  transfer  to  first  order  in  this  range.  If  the  crystals 
are  large  enough  such  that  r  >  tj,  then  the  Frossling equation,  which  was  developed  specifically  for 
spheres,  is  applicable.  To  investigate  the  influence  of  the  disk  shape  on  the  heat  transfer  rates  pre¬ 
dicted  by  the  Frossling  equation,  Wadia  ( 1 975 )  conducted  a  series  of  experiments  in  which  he  meas 
ured  the  transport  front  disks  of  various  aspect  ratios.  The  results  are  shown  in  Figure  6.  Wadia 
presented  his  results  as  the  ratio  of  the  measured  heat  transfer  coefficient  of  the  disk  edge  he  to 
the  heat  transfer  coefficient  of  the  disk  face  hf.  The  ratios  are  plotted  against  the  nondimcnsional- 
ized  values  of  r/p. 

From  his  results,  two  trends  can  be  observed.  At  low  levels  of  r/rj.  the  ratio  /te/hf  increases 
rapidly  as  the  aspect  ratio  (thickness/face  diameter)  departs  from  unity.  Wadia's  explanation  of 
this  result  was  that  the  local  shear  of  the  fluid  was  higher  near  the  edge  of  the  disk  than  on  the 
face.  For  a  fixed  aspect  ratio,  as  r  is  increased  the  he  ‘hf  ratio  decreases,  approaching  a  value  of  1 
asymptotically.  Wadia  speculated  that  as  the  face  dimension  is  larger  than  the  edge  dimension,  the 
face  of  the  disk  can  respond  to  a  larger  wave-number  range  of  the  small  eddies,  which  can  enhance 
surface  renewal  of  the  disk  boundary  layer.  It  is  at  the  larger  values  o f  r  that  turbulent  intensity 
plays  a  major  role  in  influencing  transport.  Therefore,  the  effect  of  the  turbulent  intensity  would 
be  to  enhance  hf  more  than  it  would  fie.  At  the  large  values  of>,  the  effects  of  higher  local  shear 
at  the  edges  increasing  /ie  and  the  turbulent  intensity  enhancement  of  hf  would  become  of  com¬ 
parable  magnitude,  causing  the  ratio  of  h  /hf  to  approach  unity. 

From  these  results  it  is  a  direct  step  to  determine  a  general  correlation  for  the  heat  transport 
from  a  disk  by  use  of  the  usual  dimensional  groups  of  the  modified  Frossling  equation,  but  with 
re  and  rf  as  the  characteristic  length  scales.  Wadia  determined  such  a  general  correlation.  He 
found  that  the  residts  for  all  aspect  ratios  (0.21  <  v  <  1 .0)  are  well  correlated  by  a  single  line. 

This  line  is  almost  parallel  to  the  sphere  correlation  line  but  slightly  higher  on  average. 


Wadia\s  conclusion  is  that  tor  low  aspect  disks,  the  overall  transfer  coefficients,  based  <>n  a  mean 
surface  area  equivalent  radius  |r  =  (.1  s/4tt ) 1  :  |  are  10'  higher  than  lor  spheres.  And  if  the  radius 
of  the  face  and  the  edge  thickness  are  used  as  the  characteristic  length  dimensions  for  the  lace  and 
edge,  respectively,  the  results  in  dimensionless  form  agree  well  with  those  for  the  overall  transport 
from  a  disk. 


Summary 

The  Nusselt  number  for  heat  transfer  from  suspended  particles  is  shown  in  Figure  7.  The  Nus- 
selt  number  is  calculated  using  the  properties  of  water  at  0°C,  and  a  turbulent  intensity  of  0.2  is 
assumed.  The  heat  transfer  from  spheres  and  disks  is  shown.  The  relevant  size  of  the  particle  has 
been  nondimensionali/.ed  by  the  dissipative  scale  p.  Also  shown  are  the  heat  transfer  rates  that 
would  occur  from  a  particle  at  its  terminal  velocity  under  the  influence  of  gravity. 

A  method  of  determining  the  Nusselt  number  that  provides  an  intuitively  easier  means  of  seeing 
the  relative  value  of  the  actual  heat  transfer  coefficient  is  as  follows.  Let  NuT  be  the  turbulent  Nils 
selt  number  defined  as 


hr\  h  I  v3 

Nut  =  7~  “  I  \T  * 


Let  m*  =  r/x\  where  r  =  r,  r  .  rf,  or  r  (radius  of  a  sphere).  The  heat  transfer  relationships  are  then: 
for  m*<  1/(Pr)‘  2 

Nut  =(l/w*)  +  O.I7Pr'  2 
for  1  /(Pr ) 1  2  <w*O10 

Nut  =  |(l/m*)+0.55(  Prim*)'  3  ] 
for  m*  >  1 .  with  a  low  intensity  aTw*4  3  <  1000 
Nut  =  [(l/m*)+0.70a“  OJ5(Pr/w*)1  3  ]  (3 


Figure  7.  Noiuiinwnsional  heat  transfer  correlation. 


Figure  8.  ;\<>ndimensional  heat  transfer  correlation  based  oil  a  turbulent  Nusselt 
number. 


for  m*  >  1 .  with  a  high  intensity  or,  m*4  >  1  000 

Nu,  =  |(  I  /m*)+0.70a°y25  Pr1  '|0 

where  0  =  1 .0  for  a  sphere  and  1 .1  for  a  disk.  These  Nusselt  number  relationships  are  shown  in 
Figure  8. 


NUCLEATION 
Initial  nucleation 

The  initiation  of  the  transformation  of  a  substance  from  a  less  stable  to  a  more  stable  phase  is 
called  nucleation.  When  the  temperature  of  water  is  below  the  freezing  point,  the  water  is  super¬ 
cooled  or  undercooled,  and  it  is  less  thermodynamically  stable  than  water  that  is  ice.  The  nuclea¬ 
tion  of  supercooled  water  is  the  necessary  first  step  in  the  formation  of  ice. 

Researchers  have  thought  that  frazil  ice  may  form  by  three  types  of  nucleation:  homogeneous 
nucleation.  heterogeneous  nucleation  and  secondary  nucleation.  Nucleation  in  an  absolutely  pure 
liquid  is  homogeneous  nucleation.  Nucleation  that  results  from  the  presence  of  foreign  particles  is 
heterogeneous  nucleation.  Secondary  nucleation  results  irrespective  of  its  mechanism,  only  be¬ 
cause  of  the  presence  of  ice  crystals  in  the  supercooled  liquid. 

The  importance  of  secondary  nucleation  to  frazil  ice  formation  has  long  been  recognized  (Alt- 
berg  1  036).  It  is  begging  the  question,  however,  to  use  secondary  nucleation  as  the  entire  explana¬ 
tion  tor  existence  ol  trazil  ice.  Undoubtably,  secondary  nucleation  plavs  the  major  role  in  increas¬ 
ing  the  total  numbers  of  frazil  crystals:  it  will  be  discussed  later.  The  object  heic  is  to  discuss  the 
source  of  the  original  frazil  crystals.  The  original  crystals  added  to  industrial  crystallizers  to  begin 
secondary  nucleation  are  called  seed  crystals.  Are  seed  crystals  necessary  to  start  the  formation  of 
tra/il  ice  in  natural  water  bodies?  It  so,  where  and  how  do  they  originate'.’ 

Until  fairly  recently,  researchers  didn’t  think  that  seed  crystals  were  necessary.  At  first  it  was 
thought  that  there  was  spontaneous  homogeneous  nucleation  (Barnes  1928).  It  is  well  known  today 
that  there  is  no  spontaneous  nucleation  in  pure  water  unless  the  water  is  supercooled  to  approxi¬ 
mately  -38  ('.  a  temperatme  never  found  in  any  watei  body.  Heterogeneous  nucleation  was  next 
otlered  as  the  mechanism  ( Altberg  I  636).  and  the  experiments  of  Dotsev  ( I'MX)  lent  strong 

t: 


support  to  this  ulca.  Doisey  denionstiated  that  uo  was  nucleated  In  "motes"  m  the  stipeicooled 
water,  lie  defined  motes  as  treelv  suspended  parlieles  ot  every  kind  and  mote  genciallv  as  pat  tieu- 
lar  points  ol  10  tilth  ness  in  an\  solid  boundary  that  touched  the  super  cooled  thud.  Sehaetei  (  1 9  50 ) 
tepoited  seeing  motes  m  the  tia/il  civstals  he  eiew  in  tia\s. 

I  he  chemical  and  pin  steal  characteristics  ot  the  motes  such  as  shape .  polai  Uv  and  la t  tic  mis 
match  between  lire  motes  and  ice  crystals  are  important  in  dctci  mining  their  ellccineness  in  nu¬ 
cleating  ice.  Doisey  demonsiiated  that  once  the  tcmpeiatuie  at  which  a  sample  ol  wale:  lio/e 
spontaneousK  was  determined  (or  more  accurately  .  the  temperature  at  which  the  motes  contained 
in  the  sample  caused  nucleation ).  thewatei  possessed  complete  theima!  stability  tluoiighoiii  the 
range  of  supercooling  above  that  temperatuie.  He  lound  that  samples  ol  water  could  be  super¬ 
cooled  to  within  I  t  of  their  predicted  nucleation  temperatuie  and  remain  there  lor  an  mdelmite 
tune  without  liee/mg.  However,  if  the  temperature  was  lowered  the  I  ('.  the  sample  would  liee/e. 
Any  solid  substance  can  be  charactei i/ed  by  the  minimum  supercooling  at  which  il  will  nucleate 
ice.  A  list  of  the  nucleating  temperatures  associated  with  many  substances,  organic  and  irioigamc. 
is  provided  bv  Hobbs  ( 1 074 ).  The  weakness  of  the  heterogeneous  theoiv  is  revealed  In  such  a  list . 
There  are  no  known  substances  that  will  nucleate  ice  at  the  small  levels  ot  supercooling  measured 
in  natural  water  bodies.  Piotrovieh  ( 1956)  is  credited  with  first  pointing  this  out. 

It  is  not  possible  to  state  categorically  that  a  substance  that  can  nucleate  ice  at  the  levels  of 
supercooling  measured  in  natural  water  bodies  (less  than  1  (')  does  not  exist,  only  that  none  has 
been  found.  Such  a  substance  was  not  present  in  the  water  samples  tested  by  Osteikamp  and  (lil- 
llllian  ( 1975).  which  were  taken  from  a  stream  that  was  producing  frazil.  The  samples  were  cooled 
until  they  spontaneously  froze.  The  range  of  temperature  at  freezing  was  -4.9  to  -I  5.9  (  .  The 
supercooling  of  the  stream  that  the  samples  were  taken  from  was  less  than  I  That  the  samples 
froze  over  a  range  of  temperatures  can  be  explained  by  assuming  a  random  distribution  of  motes 
between  the  samples  and  by  assuming  that  many  different  types  of  motes  were  present  in  the 
stream  water.  There  were  no  substances  contained  in  the  samples  that  could  effectively  cause  nu¬ 
cleation  at  the  supercooling  level  known  to  exist  in  the  stream.  This  is  consistent  with  the  experi¬ 
mental  results  of  Dorsey  (1948),  Margolis  ( 1 969).  Kane  (1971 ).  Ivans  ( 1 973 ),  Woltz  1 1 975 ). 

Muller  ( 1978).  and  others.  This  suggests  that  heterogeneous  nucleation  cannot  be  the  mechanism 
responsible  for  the  initial  formation  of  the  frazil  crystals. 

To  remeuy  this  weakness  in  the  theory  of  heterogeneous  nucleation.  Michel  (1971  )  has  proposed 
the  existence  of  a  thin  layer  of  very  supercooled  water  immediately  at  the  air/watei  interface.  The 
supercooling  of  this  layer  would  be  sufficient  to  cause  heterogeneous  nucleation.  somewhere  in  the 
ramie  of  4.5  C.  It  is  well  known  that  a  temperature  gradient  does  exist  at  the  .surface  of  cooling 
water  bodies  (Paulson  and  Simpson  1981 ).  However,  in  any  water  body  containing  the  slightest 
degree  of  turbulence.  It  does  not  seem  likely  that  the  surface  could  be  substantially  colder  than 
the  hulk  of  the  water. 

Observations  of  frazil  formation  indicate  that  the  first  frazil  crystals  arc  often  seen  neat  the 
wjici's  surface  (Schaefer  1950.  Arden  and  Wjglc  1972).  However,  measurements  ol  the  surface 
supercooling  level  by  thermometers  (  Arden  and  Wiglc  1 972 )  and  radiation  thermometers  (Osier- 
kamp  and  (iillillian  1975)  tailed  to  measure  inv  supercooling  level  bevtmd  a  few  tenths  of  a  degtee 
more  Ilian  that  ol  the  bulk  temperature.  Osterkamp  ( 1977)  summarized  (he  available  data  and 
concluded  that  spontaneous  heterogeneous  nucleation  in  a  thin  supercooled  suifacc  l.iyet  ol  watei 
cannot  explain  liazil  ice  nucleation.  lie  also  added  that  the  possibility  ol  heleiogeneotis  nucle.i- 
tton  existed  only  fot  very  calm  watei  surfaces  such  as  on  a  puddle,  pond  or  lake,  but  that  the  ic- 
suiting  ice  form  would  be  sheet  ice  ralliei  than  liazil  ice. 

All  the  available  data  indicate  that  spontaneous  micleation  ol  ice  is  not  possible  m  natuial  watei 
bodies  lb  a  t  are  pi  odi  icing  liazil;  diet  el  oie .  seed  civ  stills  are  uecess.iiv  .  1  be  seed  cr\  stals  ma\  come 
Itom  outside  the  water  body  or  horn  ice  alreadv  m  the  watei  body  .  There  ate  m.iuv  situations 
whcie  liazil  ice  lias  been  ohscived  in  w.ilets  in  which  ice  had  not  existed  previouslv  oi  wlicic  the 
ice  was  lar  I  loin  the  zone  poulticing  Ira /if  Oslei  kamp  (  I 7 )  pm  posed  a  mass  exchange  acioss  the 
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an  water  mteilaeo  as  the  meehaniMii  piovidmg  the  seed  ei\stals  in  these  <.jwv 
are  lee  ei  \  stals  that  have  a  \  a  nets  i*t  origins  \\  ate  I  that  01  iginaled  in  the  w  a  lei  bode  and  w  as  in 
trod  need  into  the  an  In  hnhhle  hni  sting,  splashing,  wtiidspiav .  evapoi.niuM .  ete  .  .an  t  ice/e  and 
return  in  the  I'm  in  ot  tei  ei  v  si  a  Is.  It  is  interest  me  to  note  that  otten  niininiuin  an  teinpeia  t  un's  o' 

to  -N  C  are  re  pm  ted  as  necessary  lot  the  production  o  I  ira/il  I  liese  teinpeia  lutes  cm  icsp.  aid 
to  the  minim  tun  teinpeia  tines  at  w  Inch  spontaneous  lietei  ogene.  mis  nucleatioii  c.  mid  he  e  spec  led 
in  water  pat  tides  suspended  in  air .  Ice  pat  tides  that  ot igmated  at  some  distances  1  mm  the  wale t 
hodv.  such  as  snow .  host,  ice  particles  liom  trees,  sliiuhs.  etc.,  could  he  cllcctue  seed  civstals 
Very  cold  soil  particles  and  cold  organic  materials  at  tempeiatuies  less  than  the  supeicoolmg  level 
necessary  to  cause  spontaneous  iiticleation  can  also  he  introduced  actoss  the  an  walet  mtetlace  and 
may  serve  to  nucleate  ice.  although  their  effectiveness  is  not  known. 

In  summary,  the  tm  million  of  fra/il  is  started  by  the  introduction  of  seed  crystals  into  supet- 
cooled  water.  The  mass  exchange  at  the  air/watei  interface  ptoposed  In  Ostetk.imp  ( ll,7~)  is  the 
most  probable  mechanism  In  which  the  seed  crystals  arc  introduced.  The  origin  ot  the  seed  civs- 
tals  and  the  rate  at  which  they  arc  introduced  will  depend  on  the  local  environmental  conditions. 

The  mass  exchange  mechanism  provides  a  reasonable  explanation  for  (lie  observation  ot  ice  erv  stals 
at  the  water's  surface  at  the  start  of  fra/il  formation.  This  mechanism  cannot  explain  the  existence 
ot  all  fra/il  crystals,  however.  The  number  of  ice  crystals  increases  rapidly  when  a  crystal  is  Intro¬ 
duced  into  turbulent  supercooled  water  in  which  spontaneous  nucication  is  not  possible.  This  in¬ 
crease  in  the  number  of  crystals  occurs  only  because  ot'  the  presence  ot'  the  original  seed  crystal 
secondary  nucication .  Therefore,  to  determine  the  rate  of  increase  of  fra/il  ice  crystals,  the  rate  of 
introduction  of  new'  crystals  and  the  rate  of  secondary  nucleatioii  must  be  known.  The  relative 
magnitude  of  these  rates  will  depend  on  the  local  environment;  however,  the  rate  of  secondary  tut- 
cleation  is  probably  much  greater  than  the  rate  of  introduction. 

Secondary  nucleation 

Introduction 

The  processes  that  govern  the  rates  of  secondary  nucleation  are  poorly  understood.  However, 
a  partial  modeling  of  the  kinetics  of  secondary  nucleation  is  possible  based  on  the  work  of  livans 
et  al.  ( 1974a,  b )  who  demonstrated  that  for  ice  the  production  rate  of  potential  nuclei  of  new  crys¬ 
tals  and  their  removal  from  the  parent  crystals  could  be  uncoupled.  The  most  widely  accepted 
source  of  the  potential  nuclei  is  surface  irregularities  that  are  sheared  from  the  surface  of  the  parent 
crystals  (microattrition).  Two  general  mechanisms  of  removal  of  the  nuclei  from  the  surface  of  the 
parent  crystals  have  been  suggested;  collisions  of  the  crystals  with  hard  surfaces  (including  other 
crystals)  and  fluid  shear.  If  the  rate  of  secondary  nucication  is  limited  by  the  production  rate  of 
potential  nuclei,  increases  in  the  number  of  collisions  of  an  individual  crystal  will  not  increase  the 
production  of  new  crystals.  If  the  rate  of  secondary  nucleation  is  removal-limited,  however,  the 
parent  crystals  will  produce  the  same  number  of  new  nuclei  each  collision,  independent  of  the  crys¬ 
tal’s  time  history.  Front  their  experimental  work,  livans  et  al.  ( 1974a.  b)  concluded  that  the  second¬ 
ary  nucleation  of  ice  was  limited  by  the  rate  at  which  potential  nuclei  were  removed  from  the  crys¬ 
tal  surface.  Therefore,  it  was  possible  to  determine  the  overall  nucleation  rate  ,VT .  with  two  or 
more  mechanisms  of  removal,  as  the  linear  sum  of  the  actual  nucleation  rate  attributable  to  each 
mechanism  of  removal  ( . V; ) , 

-V,  =  ,V,+,V2 +..*,■  (l>7) 

Based  on  the  work  of  (  Ion I /  and  McCabe  t  ll>72).  ( iaiabedi.iii  and  Strickland -Const able  (  1972). 
livans  et  al.  ( 1 974a.  h ),  Wolt/  (1975)  and  others,  we  will  assume  that  the  removal  of  potential  nu¬ 
clei  is  caused  solely  by  the  shear  produced  during  collisions  ot  the  parent  crystals.  The  nucleation 
rate  of  each  mechanism  ot  collision  can  he  expressed  as  the  product  of  three  functions  (Botsaris 
1970 
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tv 


•Vr  =  (F,)(F,  )(/•,) 


r>M 


where  F,  =  rate  of  energy  t  ranste  i  to  crystals  In  eolliston 

F,  =  number  of  particles  generated  pei  unit  ni  collision  envies 
F2  =  tiaetion  of  particles  surviving  to  become  nuclei 

At  this  time  the  values  ot  /•',  and  F2  must  he  determined  cm  pineal  l\ .  I  heiel"ie .  n «  -.i  i  n  p  h  t  \  mat¬ 
ters  let  f\  and  F2  he  combined  and  eq  c>x  be  rewritten  as 

•v,  VN 

where  -S’N  =  (F,  l(F2 ). 

We  expect  that  ,S'N  is  a  function  ot  all  the  parameters  that  govern  the  suitace  m.*i ph* ’!•  »c>  and  the 
crystal  growth,  including  supercooling  0 .  impurity  concciitiatioii  (  t  .  tmluilence  level  <  ,  etc.  I  he 
total  nucleation  rate  can  he  expressed  as 

\'T  =  , etc. )(/•',,  +  F,;  +/,,  .  ).  1  lull 

Tlie  next  parts  of  this  section  will  focus  oil  determining  the  rale  of  eneigv  tr.msfei  loi  each 
mechanism  of  collision.  Two  general  classes  of  crystal  collisions  can  he  identified  collisions  be¬ 
tween.  crystals  in  suspensions  (crystal-crystal  collisions)  and  collisions  between  civstals  and  exter¬ 
nal  boundaries  (crystal-boundary  collisions). 

We  will  assume  that  the  nuclei  produced  by  collisions  are  effectively  Jt  /ero  si/e.  Therefore, 
the  internal  coordinate  of  the  parent  crystal  remains  unchanged  during  a  collision .  We  also  assume 
that  large  scale  breakage  of  crystals,  which  has  not  been  observed  in  natural  watei  bodies  or  agi¬ 
tated  crystallizers,  does  not  happen  during  a  collision. 

Crystal-crystal  collisions 

We  will  assume  all  crystal-crystal  collisions  to  be  between  two  crystals  only  .  Three  types  ot 
collisions  between  crystals  can  be  identified : 

1 .  Collisions  from  the  crystals  moving  with  the  fluid.  These  collisions  aie  caused  In  spatial 
variations  in  the  fluid  motion. 

2.  Collisions  from  the  crystals  moving  relative  to  the  fluid.  These  are  collisions  caused  In 
the  crystals’  inertia.  It  was  shown  in  the  Ice  Crystal  drou  th  Rates  section  that  the 
velocity  of  the  crystal  caused  by  inertia  forces  was  always  much  less  than  the  velocity 
determined  by  the  shear  rate.  Therefore,  collisions  caused  by  inertia  will  not  be  con¬ 
sidered  further. 

3.  Collisions  caused  by  buoyancy. 

Note  that  only  collision  type  I  will  cause  collisions  between  crystals  of  similai  si/e. 

The  rate  of  energy  transfer  to  the  crystals  by  collision  can  be  determined  following  the  method 
of  Evans  et  al.  (Id74b).  F,  is  the  product  of  the  collision  energy  F(r,.  r2 )  and  the  frequency  ot 
collision  q(r,,r2)  between  crystals  of  size  r]  and  r2  .  integrated  over  the  crystal  size  distribution 
Thus 


F,  =  /  J  </{r,.r2) /(r,.  r2  )n(r,  )n(r2)drl  dr2  .  (101) 

o  o 

For  simplicity,  the  crystal  size  distribution  function,  n,  will  be  written  as  a  (unction  ot  the  cry  si al 
size  only.  We  determine  the  collision  frequency  per  unit  volume  between  crystals  of  size  /q  to  ;q  +,/tq 
and  size  r2  to  r2+Jr2  as  follows.  Let  R  =  (/q+r2)bc  the  “collision  radius."  and  '!'(/,.  r2  )  be  the 
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collision  efficiency  of  crystals  ut  si/e /,  and  r2.  The  collision  efficiency  is  defined  as  the  portion 
of  crystals  that  would  have  collided  it  the  llnid  flow  field  hadn't  been  disputed  h\  the  cisstals 
(Saffman  and  Turner  1  ‘> 5 f > ) .  Assume  that  the  coordinate  system  is  centered  on  one  panicle  and  is 
moving  with  the  steady  velocity  of  the  particle.  Let  r(r , ,  r2 )  be  the  relative  motion  between  the 
crystals  in  this  coordinate  system.  Thus 


t/O,  ,r2 )  =  n(r  ,+r2 )2  ‘(r, .  )  '\'(ry.r2 ) 


no:) 


and  by  this  definition,  r(  rx )  >  0. 

We  can  estimate  the  collision  energy  by  determining  the  energy  exchanged  between  crystals  ot 
sizes  ry  and  r2  during  a  collision.  Again,  take  one  crystal  to  be  the  origin  of  the  coordinate  sy  stem, 
say  r2 .  The  relative  velocity  between  rx  and  r2  will  be  iirx,r2).  When  the  crystals  collide  they 
will  deform  plastically  and  elastically;  their  centers  will  come  together  and  approach  a  minimum 
distance.  We  will  assume  that  the  collision  is  inelastic  enough  so  that  when  the  distance  between 
the  centers  of  the  crystals  is  a  minimum,  the  crystals  will  not  be  moving  relative  to  one  anothei . 
They  will  then  be  moving  together  at  some  resultant  velocity  r.  A  momentum  balance  is 


m( r ,  )  r(r,,  r2 )  =  l/uffi)  +m(r2  )|  r 


( 1 03 ) 


where  m(r ,  )  and  m(r2  )  are  the  masses  of  the  two  crystals  with  size  r ,  and  r2  respectively.  The 
magnitude  of  the  energy  expended  in  the  collision  will  be  equal  to  the  difference  in  the  kinetic 
energy  of  the  two  crystals  before  and  after  the  collision,  or 


h\rx.r2 )  =  ‘/a  m(rx )  r2(r,  ,r2 )  -  M  \m(ry )+  m{r2  )| 


(104) 


Substituting  from  eq  103,  we  find 


/-(/-,  ,r2 )  =  Vi  - 


m(rl)ni(r2)  2 


nt(r,  )  +  m{r2 ) 


v‘(rx.r2). 


(105) 


The  rate  of  energy  transfer  can  now  be  expressed  as 
tin  ,  >n(r ,  )m(r2 ) 

li  =  /  /  -y  (ri+r2)2  _["(7|  )77w(;2)j  T(z,.z2)r  (ry.r2)n(rl  )n{r2)dr,  dr2  (l()b) 

and  it  now  remains  for  us  to  determine  the  relative  velocity  and  collision  efficiency  applicable  to 
each  mechanism  of  collision. 

Collision  a  from  the  crystals  mining  with  the  fluid.  Lor  this  analysis  we  will  assume  that  the 
crystals  exactly  follow  the  fluid  motion,  that  is,  all  inertial  buoyancy  and  other  effects  are  ignoied. 
As  in  the  h  e  Crystal  (Irowth  Rates  section,  it  will  be  necessary  to  analyze  the  dissipative  and  iner¬ 
tial  subranges  separately. 

If  Af  <  rj.  then  we  can  estimate  the  rate  of  energy  transfer  assuming  that  both  crystals  are  in  the 
dissipative  subrange.  The  mean  square  relative  velocity  between  two  points  separated  by  a  distance 
R ( .  where  R c  <"  p  ,  is 

i  '{r,,r2 )  =  0.13  RJe/a)'  2 .  (107) 

It  R  ■  t\.  then  the  relative  velocity  cannot  be  estimated  as  neatly  as  above.  It  may  be  that  /•,  o 
r2  or  both  arc  greater  than  r/.  We  will  assume  that  if  R  >  rj.  then  the  crystals  will  move  with  a  tela 
tive  velocity  appropriate  to  the  inertial  subrange,  regardless  of  each  individual  crystal's  size.  The 
mean  square  velocity  difference  between  points  separated  by  a  distance  R  where  R  •  rj  is 


v(rl,r1)  =  2J(eRi:)'  \ 


(ll)M 
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The  collision  efficiency  ,r2 )  is  difficult  to  estimate.  Levicli  (  I  9r>2 )  discusses  the  pioblem 
of  determining  the  collision  efficiency  of  suspended  particles.  Complications  arise  because  neither 
particle  is  stationary,  because  turbulence  and  the  boundary  layers  around  the  particles  modilv  the 
fluid  stream  lines,  because  inertial  forces  on  the  particles  must  be  accounted  tor.  etc.  I’ruppaclier 
and  Klett  (1978)  discuss  the  uncertainties  of  calculating  the  collision  efficiencies  ot  water  drops  in 
air.  From  experimental  evidence,  Saffman  and  Turner  (1956)  assumed  that  the  collision  efficiency 
was  unity  for  collisions  between  particles  of  nearly  equal  si/e  when  AIL  <  r/.  llowevet.  the  discus¬ 
sions  of  Pruppacher  and  Klett  (1978)  indicate  that  this  may  not  always  be  true.  In  general,  the 
collision  efficiency  goes  to  zero  as  r,/r2  goes  to  zero.  As  rx  approaches  r 2 .  the  collision  efficiency 
rapidly  rises  to  a  maximum  which  may  equal  unity,  and  then  may  increase  or  decrease  as  rx  becomes 
exactly  equal  to  r2 .  In  short,  estimation  of  collision  efficiency  must  be  an  uncertain  undertaking. 
Therefore,  we  will  arbitrarily  set  the  collision  efficiency  equal  to  unity,  but  not  without  reserva¬ 
tions.  When  r2  >  r , ,  the  energy  of  the  collision  will  be  proportional  to 


m(r,)m(r2) 

m(rl)+m(r2) 


m(r ,) 


(109) 


and  so  will  be  very  small,  as  m(rx )  =  rx  .  Therefore,  the  error  introduced  by  u  t  having  'I'lrq .  r 2 ) 
go  to  zero  as  r , /r2  goes  to  zero  should  be  small. 

Substituting  the  above  expressions  into  the  rate  of  energy  transfer,  we  see  that  for  R  <  r? 

Et  =  (0.003 )(e/u)3  2  f  j~  (r,  +r2 )5  |  "(r, )  >‘(r2  )drxdr2  (II 0) 

and  forf?c  >  r? 

=  (30.92)(e)  / 1  (r,  +r2 )3  "{r'  » "{r>  )dr'J^  '  ( 1 1  1  > 

Collisions  caused  by  buoyancy.  When  the  density  of  the  crystals  is  not  identical  to  that  of  the 
fluid,  the  crystals  will  move  relative  to  the  fluid  under  the  influence  of  gravity.  This  motion  will 
be  in  addition  to  all  other  motions  of  the  crystals.  In  the  Heat  Transfer  from  Ice  Crystals  Sus¬ 
pended  in  Turbulent  Water  section,  the  relative  velocity  induced  by  gravity  was  determined.  Using 
the  intermediate  law  for  calculating  the  drag  coefficient,  we  can  find  the  terminal  rise  velocity  of  a 
crystal  to  be 


Ut  =  (0.16^'°  7,5/iz0*42*>rl  14 


(112) 


and  the  relative  velocity  v(r , ,  r2 )  is  therefore 

v(ti,r2)  =  (0.16g’o'7,5/po'428)(|r  i1'14  -rj,l4|). 
Again  collision  efficiency  will  be  equal  to  unity.  Thus 
£,=(0.16 g,OJ'5lvOA2H)\rtl2)  1 1  (r,-hr2) 


(113) 


o  o 


m(r,  )m(r2 ) 
\m(rl  )+m(r2)| 


(|r,ll4-r2,l4|)1  n(rt )  n(rx)  n(r2)  drx  dr2. 


(114) 


Crystal-boundary  collisions 

The  second  general  class  of  collisions  are  those  between  the  crystals  and  the  external  boundaries. 
To  determine  the  frequency  of  crystal-boundary  collisions  requires  knowledge  of  the  si/e  and  shape 
of  the  water  body  of  interest.  As  this  report  is  intended  to  be  as  general  as  possible,  tins  knowledge 
cannot  be  assumed  ,  however,  this  section  is  presented  in  the  interest  of  completeness. 

The  experiments  of  Evans  et  al.  ( 1974b)  demonstrated  that  collisions  between  ice  crystals  and 
the  walls,  baffles  and  impeller  of  the  crystallizer  were  a  significant  cause  of  secondary  nucleation. 
They  assumed  that  any  crystal  moving  with  the  bulk  How  could  potentially  collide  with  the  im¬ 
peller.  and  any  crystal  closer  than  an  eddy  size  away  from  the  wall  could  collide  with  the  wall. 

They  were  unable  to  determine  if  either  collision  mechanism  dominated,  but  they  did  find  that 
coating  the  metal  surfaces  of  the  impeller  and  crystallizer  with  a  soft  material  substantially  re¬ 
duced  the  nucleation  rate. 

In  natural  water  bodies,  collisions  with  external  boundaries  are  probably  not  a  significant  cause 
oi  secondary  nucleation  except  in  some  circumstances.  First  of  all.  the  ratio  of  surface  area  to  vol¬ 
ume  ot  most  natural  water  bodies  is  much  smaller  than  that  of  a  crystallizer.  Also,  the  boundaries 
of  natural  water  bodies  tend  to  be  very  rough  compared  to  crystallizers.  Crystals  colliding  with 
these  rough  boundaries  tend  to  stick  and  remain  at  the  boundary,  which  results  in  the  buildup  of 
anchor  ice.  Therefore,  only  in  locations  such  as  shallow,  rocky  rapids  could  crystal-boundary  col¬ 
lisions  be  important. 

The  rate  of  energy  transfer  during  the  collision  of  a  crystal  with  a  boundary  may  be  crudely 
estimated  in  the  following  manner.  This  analysis  does  not  account  for  the  effects  of  a  boundary 
layer,  inertia,  etc.  Let  (!m  be  the  maximum  eddy  size  that  can  bring  crystals  into  contact  with  the 
boundary.  Assume  that  the  number  density  of  crystals  is  uniform  throughout  the  region.  The 
probability  that  a  crystal  of  size  /  exists  in  the  region  within  a  distance  close  enough  to  collide  with 
a  boundary  is 


51  jr1  I  "(r)  Jr 


where  /?h  is  the  hydraulic  radius,  or  the  ratio  of  volume  to  surface  area.  The  relative  velocity  of 
the  crystal  and  the  boundary  can  be  estimated  as 

^rb  *  ej1'3  (IK 

and  the  energy  at  collision 

/trb  =  \i  m(r)  i;2h  =  Vi  m(r)  (e  CJ23.  (Ill 

Now  the  rate  of  collisions  will  be  proportional  to 
ie  C  V  -1  <i  r 

«rb  -f  I  «<rf  c/r  (Ilf 


and  thus 


e^n_  f 

2  Rh  J 


m(r)  n(r)  dr. 


Summary 

In  this  section  (he  kinetics  of  secondary  nucleation  has  been  partially  modeled.  Starting  from 
a  theoretical  formulation  of  the  rate  of  secondary  nucleation,  we  find 


,Vt  =  (^«,+£,2+/;-,3+...)XN(fl.e.Cr.ctc.). 


Tlie  rate  of  energy  transfer  Ftj  has  been  determined  for  several  different  mechanisms  of  collision. 
The  function  SN  is  necessarily  empirical  at  this  time  because  we  know  of  no  way  to  calculate  its 
value  theoretically.  i'N  is  the  product  of  two  functions:  h\ ,  the  number  of  particles  produced  per 
unit  of  collision  energy  and  F2 ,  the  number  of  particles  surviving  to  become  crystals.  We  expect 
that  F,  should  be  a  linear  function  of  the  collision  energy  and  have  relatively  little  dependence  on 
the  supercooling  or  turbulence  levels.  F2  should  largely  depend  on  the  supercooling,  in  accordance 
with  the  survival  theory  (Garabedian  and  Strickland-Constable  1974).  .S'N  should  largely  depend 
on  the  supercooling  and  perhaps  depend  less  on  the  level  of  turbulence,  e.  The  experimental  re¬ 
sults  of  F.vans  et  al.  ( 1974a,  b)  are  difficult  to  analyze  conclusively.  However,  they  reported 
that  the  dependence  of  SNon  e  is  small. 

The  terms  in  eq  1  10,  1 1  1  and  1 14  for  l:\  are  difficult  to  compare  with  each  other.  To  facilitate 
a  comparison,  let  r  be  the  average  crystal  size,  specified  in  any  convenient  manner.  Then  let  r*  be 
the  nondimensional  crystal  siz.e,  defined  as 


We  will  assume  that  the  crystal  size  distribution  can  be  expressed  as 

n{r)  =  /V(0)  T(r/r)  (121) 

where  T(r/r)  defines  the  form  of  the  crystal  size  distribution  and  is  constant  with  time;  /V(0)  is  the 
number  of  crystals  of  a  size  r  =  0.  When  the  form  of  the  size  distribution  is  fixed,  specifying  the 
number  of  crystals  of  any  size  determines  the  number  of  all  other  sizes.  Expressing  the  size  dis¬ 
tribution  in  this  way  is  a  simplification  introduced  purely  as  a  device  to  facilitate  comparison  of 
the  mechanisms  of  collision.  It  may  not  be  possible  to  do  this  in  actual  practice,  as  the  form  of 
the  size  distribution  may  vary  with  time. 

For  example,  introducing  eq  120  and  121  into  1 10  and  noting  m(r)  =  p,Kvr3 ,  we  see  that 

=(0.003  )(e/v)3  2PiKvP°N(0)2f  [  (rf  +  rff  T(rt)T(rl)Jr,  dr2.  (122) 

o  o  12 

Inside  the  integral  are  only  nondimensional  terms  and  the  value  of  the  integral  will  be  a  constant 
with  time.  Combining  all  constants,  we  can  now  write  eq  1 22 

£,  =  ai(e/u)3,1r8  0  (123) 

and  in  the  same  manner  we  can  write  eq  1 1 1  and  1  14  for  Rc  >  r? 

=a2  er6'°  +  a}rHA2.  (124) 

The  first  term  on  the  right  side  of  each  expression  is  the  rate  of  energy  transfer  due  to  crystal- 
crystal  collisions.  The  second  term  in  eq  124  is  that  due  to  crystal-crystal  collisions  caused  by 
gravity.  The  rate  due  to  crystal-boundary  collisions  has  not  been  included  here,  as  its  determina¬ 
tion  requires  knowledge  of  the  particular  water  body.  The  magnitudes  of  Et  are  nonlinear  functions 
of  the  form  and  magnitude  of  the  crystal  size  distributions.  It  is  easy  to  see  that  small  changes  in 
the  crystal  size  distribution  could  cause  very  large  changes  in  the  rate  of  secondary  nucleation. 

The  sudden  appearance  of  vast  numbers  of  frazil  crystals,  observed  both  in  natural  water  bodies 
and  experimentally,  may  be  a  result  of  this  large  nonlinearity  in  the  rate  of  secondary  nucleation. 


FRAZIL  ICE  DYNAMICS 


Introduction 

In  this  report  the  mathematics  of  the  formation  of  fra/il  ice.  the  basic  equations  and  expressions 
for  the  important  parameters,  were  developed.  Much  work  remains  to  be  done.  The  author  hopes 
that  the  synthesis  contained  in  this  report  will  provide  insights  into  the  phenomenon  of  fra/il  and 
provide  a  foundation  for  further  work. 

This  section  presents  a  simple  but  practical  example  of  the  crystal  number  continuity  equation 
and  the  heat  balance. 

Steady-state  crystal  number  continuity  equation  and 
heat  balance  in  a  MSMPR  crystallizer 

Equation  9  can  be  specifically  and  practically  applied  to  a  mixed  suspension,  mixed  product  re¬ 
moval  (MSMPR)  crystallizer  operating  at  steady  state.  This  application  results  in  the  most  widely 
known  and  used  analytical  solution  to  eq  9.  Because  of  the  generality  of  eq  9,  it  is  not  in  the  best 
torm  for  this  particular  application.  A  very  usable  form  can  be  achieved  by  integrating  this  equa¬ 
tion  over  the  macroscopic  volume  ('(/)  in  the  external  coordinate  subregion  that  corresponds  to  the 
volume  of  the  MSMPR  crystallizer.  To  carry  out  this  integration,  the  volume  of  the  crystallizer  must 
be  well  mixed.  That  is,  in  any  arbitrary  small  element  of  its  volume,  a  full  and  representative  crystal 
population  distribution  that  is  dependent  only  on  the  internal  coordinate  must  exist.  Within  this 
volume,  the  population  density  function  n(R.  r)  can  be  written  n(r,  t )  as  it  has  no  dependence  on 
the  spatial  coordinates.  n(r)  can  now  be  interpreted  as  a  size  distribution  function.  In  addition,  O', 

B  and  D  are  also  required  to  be  independent  of  the  spatial  coordinates  within  this  volume.  There¬ 
fore  multiplying  eq  9  by  r/Kand  integrating  over  V(r),  we  obtain 


/ 


~t  +  V  (Teti)  +  j-r(Gn)+D-B  dV  =  0. 


(125) 


Every  term  can  be  removed  from  the  integral  except  for  the  second,  which  results  in  an  integral 
over  the  volume  of  the  spatial  divergence  of  the  population  flux.  This  term  can  be  converted  into 
a  surface  integral  of  the  population  flux  flowing  through  the  surface  of  the  volume  V(t ),  which  in 
principle  may  not  be  stationary.  This  results  in 


fv(?en)dV=  (126) 

V  K 

where  QK  is  the  flow  rate  and  nK  the  population  density  of  the  Kth  input  and  output  stream  into 
V(t).  Carrying  out  the  integration  of  eq  125  and  dividing  by  L(r),  we  find 


+  —~(Gn)  +  D-B+ti 
dr  dr 


c/logK(r)  _  y> 
dt  V 


(127) 


This  is  the  number  continuity  equation  integrated  over  a  macroscopic  external  volume.  It  is  aver¬ 
aged  in  the  external  phase  space  and  distributed  in  the  internal  phase  space. 

In  the  same  manner  the  heat  balance  equation  (eq  16b)  can  be  integrated  over  the  macroscopic 
volume  of  the  crystallizer.  This  results  in 

dO  _  I,  dPj  _  y.  Qk° K  ,  Qk  PiK  Q* 

dr  Cpf  dr  v  '  ^  VCpf  CPf  ll-*' 

where  pjK  is  the  mass  density  of  ice  in  the  A'th  input  or  output  stream  into  T(r). 


The  following  assumptions  regarding  the  operation  of  the  MSMI’R  ciy^lulli/ei  can  now  he  mad 

1 .  The  volume  of  the  crystallizer  does  not  change  with  rune.  I  Ins  requires  that  the  volume 
of  crystals  he  small  compared  to  the  total  volume  (a  thin  suspension). 

2.  The  remaining  time-dependent  terms  can  be  d topped  as  the  crystallizer  operates  at  steady 
state. 

3.  The  birtlr  and  death  functions  can  be  set  to  zeio  for  all  sizes  of  crystals. 

4.  The  crystallizer  can  be  assumed  to  operate  with  clear  or  unseeded  feed  streams.  Also, 
the  size  distribution  of  crystals  in  the  outflow  is  exactly  identical  to  the  distribution  m 
the  crystallizer.  The  ratio  l'(/)/l(3K  can  then  be  interpreted  as  the  flow  -through  time.  r. 

5.  McCabe's  law  is  assumed  to  hold.  This  law.  often  applied  under  the  conditions  lound 
in  industrial  crystallizers,  requires  that  the  growth  rate  not  be  a  function  of  the  crystal 
size. 

(>.  The  population  density  of  newly  nucleated  crystals  will  be  defined  as  //"  and  will  be  as¬ 
sumed  constant.  We  assume  the  size  of  the  nucleated  crystals  to  he  vanishingly  close  to 
zero,  //"  is  the  lower  boundary  condition  for  the  crystal  size  distribution.  The  relation¬ 
ship  between  rr"  and  the  secondary  nucleation  rate  A',  is 

»"  -  (.V ,/(/’)  (124 

liquations  127  and  128  can  now  be  written 


.  3  n  n 
'  dr  +  t 


=  0 


(130 


(131 


Separating  variables  and  integrating  eq  I  24.  we  obtain 

f‘((Jn/n)  =  -  f  (Jr/Gr)  (132 

n°  6 

which  gives 

//(/•)  =  n"  expi-r/Gr)  (133 

and,  from  eq  131 .  assuming  0jn  =  0  t  gives  us 

Pj  =  -  ((?*//.  )t.  (134 

Crystal  size  distributions  of  MSMPR  crystallizers  operating  at  steady  slate  are  often  very  close 
to  the  size  distributions  of  eq  133.  This  provides  a  convenient  graphic  method  of  determining  the 
growth  rates  and  the  density  of  nucleated  crystals  by  plotting  the  log  of  measured  population 
density  against  the  crystal  size  r.  The  graph  is  a  straight  line  with  intercept  of  log//"  and  slope  ol 
-(fi  r)  1  .  Additional  information  can  be  obtained  about  the  growth  rates,  the  mass  of  crystals  or 
other  aspects  of  the  crystallizer  operation,  lor  example,  noting  that 

Pi  =  P,Av  f  r' 

0 


substituting  the  distribution  for  N(r)  from  eq  \^^  and  solving  tor  (>’.  wo  get 


3| 


O 


Run  216 


\(lll"KJ.  p,  T  I 


t I  •'  3 ) 


Increasing  the  rate  at  which  heat  is  removed  or 
decreasing  the  residence  time  will  increase  the 
average  growth  rale  of  the  crystals.  It  is  in  this 
manner  that  the  operations  of  such  crvstalli/ers 
are  studied  to  determine  their  operating  charac¬ 
teristics. 

As  was  seen  in  the  Ice  Crystal  Growth  Rates 
section,  the  growth  rate  is  not  constant  with 
crystal  size.  This  will  cause  the  crystal  size  dis¬ 
tribution  to  deviate  from  the  straight  line  of  eq 
133.  If  the  growth  rate  varies  with  crystal  size, 
eq  1  is  written 


d(Gn)  +  n  _  Q 
dr  t 

This  equation  can  be  solved  as  the  difference  in 
the  product  of  the  growth  rate  and  the  distribu¬ 
tion  at  any  two  sizes 


Figure  9.  Frazil  crystal  size  distribution  and  r 

growth  rate  (after  Smith  and  Sarojim  1979).  (G n)r t  -  (G n)r2  =  -(1/r)  j  n{r)dr. 

'2 

This  expression  allows  the  variable  growth  rate  to  be  determined  from  the  measured  crystal  size 
distribution  if  the  growth  rate  at  any  size  is  known.  An  example  of  such  data  is  shown  in  Figure  9. 


SUMMARY  AND  FUTURE  WORK 
Overview 

The  equations  describing  the  dynamic  interaction  of  the  crystal  distribution  and  the  heat  balance 
of  the  water  were  developed  in  the  Basic  Fquations  section.  Randolph  and  Larson  ( 1971)  used  the 
term  “information”  to  describe  the  crystal  distribution  and  the  heat  balance.  “Frazil  ice  dynamics” 
describes  the  unique  internal  information  feedback  loop  that  regulates  the  crystal  size  distribution. 
The  mechanism  for  this  information  feedback  is  through  the  level  of  supercooling  resulting  from 
the  balance  of  the  latent  heat  released  by  the  production  of  ice  and  the  total  heat  loss  from  the 
water.  The  supercooling  level  of  the  water  ultimately  determines  the  rates  of  secondary  nucleation 
and  crystal  growth,  which  in  turn  determine  the  dynamic  crystal  distribution  at  any  time.  The 
level  of  turbulence  controls  the  rate  at  which  the  feedback  loop  can  operate.  This  is  shown  graph¬ 
ically  in  Figure  10. 

The  two  basic  environmental  parameters  are  the  rate  of  heat  loss  from  the  water  and  the  level  of 
turbulence  of  the  water.  The  rate  of  heat  loss  can  be  thought  of  as  a  disturbance  that  manifests  in 
the  supercooling  of  the  water.  The  crystal  distribution  develops  in  response  to  this  disturbance  and 
its  ultimate  effect  is  to  eliminate  the  disturbance  and  return  the  temperature  of  the  water  to  the 
freezing  point.  The  level  of  turbulence  controls  the  rate  at  which  the  crystal  distribution  can  re¬ 
spond  to  the  disturbance.  The  level  of  supercooling  is  determined  bv  both  the  rate  of  cooling  and 
the  level  of  turbulence.  Increasing  the  level  of  turbulence  will  tend  to  decrease  the  level  of  super¬ 
cooling.  This  may  have  profound  implications  for  engineers  who  design  water  intakes  that  must 
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9  Supercooled  Temperatjre 
A-  -  Total  Surface  Area 
Available  for  Growth 

Q*  -Heat  loss 
6 -Level  of  Turbulence 
'energy/mass  of  fluid) 


G- Linear  Growth  Rate 
h-Heat  Transfer  Coefficient 
Nr- Secondary  Nucleaiion  Rate 
f(n)  -  Size  Distribution  Dependence 
of  Secondary  Nudeation 


Figure  10.  Frazil  ice  dynamics. 


operate  in  winter.  This  interaction  between  the  rate  of  heat  loss,  turbulence  and  supercooling  was 
suggested  in  the  experiments  of  Carstens  ( 1966). 

Basic  equations 

The  basic  equations  describing  the  crystal  number  continuity  equation  and  the  heat  balance  of 
the  water  were  developed  in  the  Basic  Equations  section.  Simple  versions  of  the  heat  balance  have 
been  applied  to  the  formation  of  frazil  ice  in  natural  water  bodies  by  previous  researchers  (Oster- 
kamp  1978).  The  number  continuity  equation  has  not  been  applied,  and  it  is  unclear  at  this  point 
what  the  resulting  distribution  will  be.  The  size  distribution  of  frazil  crystals  Iras  not  been  measured 
for  any  water  body  under  any  conditions.  These  are  the  first  vital  field  data  that  should  be  collected. 
The  hydraulic  and  meteorologic  conditions  under  which  the  size  distribution  is  measured  should  be 
carefully  determined.  Parallel  to  this  effort,  the  solution  of  the  basic  equations  should  be  pursued. 
Currently,  the  solution  of  these  equations  cannot  be  considered  routine.  The  numerical  schemes 
necessary  for  their  solution  must  be  developed  and  tested. 

The  importance  to  the  crystal  distribution  of  such  physical  processes  as  the  formation  of  anchor 
ice  and  flocculation  of  the  frazil  disks  is  not  yet  known.  These  processes  must  be  investigated  and 
their  influences  assessed. 

Growth  rates 

Despite  much  effort,  the  intrinsic  kinetics  of  ice  crystal  growth  has  not  been  completely  deter¬ 
mined.  This  is  the  one  major  potential  problem  in  determining  the  growth  rate  of  ice.  Determina¬ 
tion  of  the  crystalline  kinetics  requires  carefully  controlled  experiments  in  which  a  water  tempera¬ 
ture  can  be  accurately  maintained  and  measured.  For  growth  controlled  by  heat  transfer,  the 
growth  rates  can  be  estimated  fairly  well.  However,  the  transfer  rates  in  the  inertial  subrange  have 
been  determined  based  on  the  Frossling  equation,  an  experimentally  based  relationship.  The  de¬ 
velopment  of  a  rational  theory  of  transfer  rates  in  this  region  is  needed. 

Nucleation  rates 

The  nucleation  rates  of  frazil  ice  crystals  are  probably  the  largest  unknown  at  this  time.  The 
existence  of  secondary  nucleation  of  frazil  ice  in  natural  water  bodies  has  been  conjectured  (Ostei- 
kamp  1977)  but  not  yet  demonstrated.  The  development  of  a  theoretical  basis  for  the  determina¬ 
tion  of  the  number  of  nuclei  produced  per  unit  collision  energy  is  needed  as  is  a  rigorous  testing 
of  the  survival  theory  and  experimental  confirmation  of  the  collision  rates. 
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